


Institutional Archive of the Naval Postgraduate School 


Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations l. Thesis and Dissertation Collection, all items 


1991-09 


A comparison of three numerical methods for 
updating regressions. 


Raptis, Grigorios J. 


Monterey, California. Naval Postgraduate School 


http://hdl.handle.net/10945/26419 


Downloaded from NPS Archive: Calhoun 


| Calhoun is the Naval Postgraduate School's public access digital repository for 
D U DLEY research materials and institutional publications created by the NPS community. 
3 Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS's first 
KNOX appointed — and published — scholarly author. 


LIBRARY Dudley Knox Library / Naval Postgraduate School 
411 Dyer Road / 1 University Circle 
Monterey, California USA 93943 





http://www.nps.edu/library 









































* 
he a A p AUG 
r D & à P T 
"^ ik AP a " e 8 S 
‘ e *. 1 LN ia ] 3 
A T 3 TA T 





E 
SS SF 


Se 


t aa tyit 
; - LAM. es ED DE A d 

" 

e 


"d 
m m S 
i Y x 
x CU A A MAE URS eL Ww Ai C CTI 


1 De CON CCP Sn LIT m. A ORAL A 
x E ' [17 E Abb ane L v 
J T JA OO BCE le otha” race 
E $ Aah ?, LEPE TS 4 Pa E t wie k ` 
he i pe% % 3 Err 3 
* * " 


ws 
b fe ^v USA T€, 
SN XX 2: ae oe 
LE M $" 5» ben." Y * 9393/4. "oe ae is «n 
ST Pappa A OLY. AA LITI We D er YN a Mite s IR SPI LL Ad ELE X CS E hà. d 
NYXM L] AL CEPS J^ pe SHE A e RC aV AG EERE CLIN ICTS yo PA v Sw 
UE ALIUM VS e LV w EE ATEM CIC earl SOO PO Gone ert TNT Bap a ea ow fr v eee 
a LEMP ys uit «UE XS ate RL Pre SERA 1 ry AI X" LLL 
= 1 POO aU We CODE e DOE [] COO eS a LZ LI 
4 | a Penh ti ee MENSES Wm AMAL LL EY Sere Ala MOLAS KF uh d Es i 
` "a CES APP" rrr A E Er SES 4. ade ^ Le EU Pda RA X ops ?] LA ECC 
f j " r *. TUE IET 1s De WDA gia N id ETT M Ak ee oy Am Nd 
p " tW ry P A ees E EN E LU pA ALIUS yt t E "a ALTUS OE ee so DAC: 
s r E Ph. [1 » * PR4aee Pei Sd o rn Ag [TR DP TS ae 
7 ANN Y PLOTS E Ty. RETA CL Ca s e dris 
Ps LTS es Let ee eT B VES PRL yy. Toys P 
z AP, LEN IN 2% 4 ey rs SACRO hg bp e RCRA CATE NN 
fy p e , 
1 ty 2A X X C or DC OO Y CREE CY aa vas A Ye VL ts het 
Zu A TE aR GL UTET XE ud USED PRP eH a weatan ir 
D AE "255 95V yu aula NX : hA 
E Lp [ eer) 


cay L] E " 
A : LS jN p ARA a tp r3 RE ik Ler n 
LI 3 TETAS Anen A rI we ^e do Eon NT no ee my KE Uu 
pw). oe pe% obe LEY Cu Po r cm PW rm fo 
p dL T -^.4 j ere P 


E 
t 
ue ve Pare d'a el 
x ERR aT AN S ex ^ ULT RC Ms USE past Ld 
aide li de tk te + a he y d Ula ML d v. 
i . MUN Pea, c9 ORAE EAD CM 5 ' SERENO 
^ TE Y ] 9? 4 LN LU Beet Fo at, nd f a o 
Sn oh d DEC E ET i A ot he Im n T v 
M | t J ei ia AA A EA ERR E Edo XC na 
m LAM ELSE Mean eye n» eh Ais LL i areta h 
7 EEE. pe AERLE N ma i A AZ 2 M E ri LL RECUPERA TM AC e 
ey SLE A E TE KE łe Le At Mel ENT fea CES aby me any Sit Pere OS® URAS PAR CURT MEC ALI 
à T ANIMO" " M. i ME Li 2 OME Ah oy © palate Hh bo i UAR C LRL ES LI SS LS bre ya) 
3 a : E SORT VAT E e AR. RN LCD DEO P oh ind a 
URL CC EC! SG LAREDO b M Lu LIE RC Ra DS 1. PERMET 
LL E X ® tre’ns ERR cn So E ees Yo he k : LE XA "el IE, Ld A! ERR ATY 
^" ys oa eee Te i 2 6 rt e“ BA LIENS ALPE Oe Orr wage Apel OC Por Dr 2 adr be TaT bie as [Lt 
" ' ie UC Ge ee ees Sr Jta OS mad US PO POU A ey UM aves cg ed ate hd Pt eh ARTY 
;9. y SO ary ak, ED A tak LU nie AG SD Ie Lae to a a EARN RR eae UR C. M 
CLE LO CN ae PY : ACIER LO V AL chad oo dee hod he LT ee tee RV. SS 
- " hee © Tae. eee rs Dh shine eM CP Yu v» 
* P" LL T LOE E O^ Ye dE LER: E 
I F* " b's 1 [P y $ Fd 
Ce 1 LOHN i | 


" 
hd] 


4 & Sead ap 
LN awh d] A RUD L dns | w 
^s d ` , 






N 4 QN E Tw dE om 
IT. d s ete wt EE TY LII LPs "uM A LIE e 
a WEIT DRESD ue bnc or «Ex at TGAN er aE INES 
; seft yee ies a y pe Ce ps PNrM rte 
PRÉ v. EE ¥e PP pt. X "T CRX $44 54 ma ve Te "m ET 
L| ATE CLP TO Piste ER Mo Pah rs Tec ed as reri as ges ti: - wn 
$ y. 

EE E £X" y : X f LER REL i XS HERD ave AM oC re ia zar Ar. SX DAT ttt rris 19 A 

M a a laa 1e MEL VEO KA AAR TRN W& Maas M [2 AU xs: LR SORTA 
E "Xy Pur" > me i dy TE emanate PAR re ys Lm d RS a 

M 9 uS * no XX d i & tr. t lh. tial, Thao D Bale M Ag. "d T un En 
N ) 18 Fay a. eS OTERAO NY ‘ea AC 

i Fh ET i ae | d x oe Arle n mA X 
DA ECL teh Me TOMAS sd A afta 

gt paw Say 3 x 


r n oe 3 (^ v. Aa Wrak MEN pret (A S2 
2 ak ee} he et TS rf S LEY et eai 
D DAN PLA hero P TAX Ls i " 
$È Ys Mes. vw 
































NU As A nem 
eb dido ELS E by? à i tae 
D E hg X m ae un X ded ^ A VR (EY Nd VU 
o QE A wh y E iM "7 a dae CS PRA Ke ORAT ACT. D [T3 4 Ms. 
VN y A rx tet Ata LM x CERE > ei POE EUN AE E MCA Ei ARE RA 
SPAR: eT Eee OE ag w art ooh AG ANM ne X e AI ra mes 
de vbt: PISA STD RU Petes Got Ny ere T Ses Lu y ra 
ty RU Niet i 0: Ems? kaos, DEL Non EN 3 TAA PANA RD Loa fa 77 
P e 79 i HN EA ARES "2 »* eri PAAR t- AADI r 3 ES i. ree P p » CODE DT 
Pe ee re B «i P am" G P X a - NARAS 
" Wia a” d hA ~ E LS; ket. Bade: QU dá Da C Ded a nA I " halve aie A h r * : ei 
22 PL rure M Mi » 


tL EU 
oh eee S XEM 1 
UT rn 2x HARE 

goes À 


bic y 

TM. pop pt aie, A tut b EN: of 
ert. ae e Mrs we d WX" I 4 La C 
? > * ¢ Ret DU Tar R eres ae n H a7 RARER ARE ARGU ACCT H rwn asi 
AM ID ee Ti LY Ih ERE are A rrr MAPS a "t T DUE S E 
; j n* HAE TA Y le - rrr A a cl 

e ys OD tae hes " i sl + kéo P vt "d ne: TuS AK xu BUE. e c4 A ttf ir EM ICT DECR ET d 

y Lu p de bo he) thd at oe PX. HS 

Li 4. n LR 


=~ 

1 (3 
1 ETIN 2) Ei rr MP gen ey 
Vae vr EON An ER CA A Pee 
SM Mee IE des A M CREE 
A s +; 05 3E arl Mx - een menu H d AS T et Yet are RA LA a OA Y 

9.7 TRS a EG VECES Tails sare P dx Mt A OE 
Z TOS MOS LS ies Ty EE OPE OS «e D CERE CE era: * 

D M A 
E 





= 















































A SOR tech Car LX Ra EX 
y x d P CE ir M En EE RS 
LET 1 ea yp a tafa NW P em e Rte UE NERO Tt ats Ig Bathe Pour * 
| PA Wi x SANI AU FATA PRAET! E DECEM EE DRE s BEIM od Ri 
PEE "ese A Ye Pay P Dt Chl ADEL RS 
E GL nn eee a Bea ERU der A 
CARRE ers A OI d. eit? Brak 
«ey riu? 7 V u^ C) Q Pira 
f INN E A SE A hey o OR 98 LA a 
^ F Qe d" "A r n ^) 4 Se bd 
A PGE RIO ARE p SE QU RAZAR NL ju Ne 
poe. abc oU DOO VAY E] Ae 
i SR Cit tise ; : EN. (S 
a AAT LAE 
4 EAS i H M 
LI fi. we ^ s 
YY a ls i 
+ yy Vie sad ub i 
aie Pun Pd E C 


h ^ 
DOCE E? 
MEC EC QUOKA DORIA 
D [LES uu. de f y ru . E ef », 
DUCERE Cun A12 E CS 
n2 72 x NC P. H 
ll 3 










qu y e 
"IMS 


at E 


uA 


a 


n 


e T 
LONE TEL 
E 








































































2 E RA E E , ia LL 
E nero eR y EAE ER 
5 mne ve DUC TN DONC ee ee te oe 
x 3 COS A 3 eA HE KI ere mt: 
E à Ma VAS P Re HE 
mr n [QEON TEE 5 
, N A AE ALAS AX SC PU ea oce Nok ety 
I MEM DELE: M E rl, EY PULO CC t NICK, 5 
^ M E fa 4 LL " MN MC ar) ) uUa n Pert L4 e : LAE 
: A ae) " O i NUS ELE, s A ad a xu A Feet RO DU DOS VETE OU CES 
g L E H LA] * Jis BLATT P P »*Ad' e á P, 
n I l ^ F me a A ` ND "wn ns. " TER UU Ne wx 3723. C$ TUO an de E 
g E » D U ua, TUS , ny Uu 3 UTE D o 
A è f : Jd vy t E' n ‘at? Fe , . DEN EAE LEE DA R% a Se rM Doreen 
" ' ID id dor ae LE DON TE E SYN ue LE On » 
r 5 , a. s MP A pase Ba Tre, gy X MT ) UE SER x dtes ma EEA 
F G A ae Pak EE UP : tv) EY Jr DE i (A reo) de NAE EU CE Rey, 
A ^ 9 Uo ONE, gat E H F KRAE E OE E Ww aes MAC j 
| d^ .e i G 1 "rm ACE M DEE TAL i 
29523 4 44 P* " i Y PE. Pha SS ae fs SARE 4 ; are, if 
LI TL * . c 
a’, i EUN EL "Lu HR : TU tats XUL y ae AN i 
o )* Ag ou ayaa vs Ley ULL TNCS 
A E TE bern pA g F A Pr y^ A6 BA, Rot eed 
C ? mm n Yn E " Y Ag ev d uta M Hg Pet Prid eut 
3 2 E "n Py st CS Mn ELM s ed Ws 3a fondo v uj OG f Pe eee ts be 
he En fe & P Ó 733234 Ha AP b vi I VIAL d SR uus Mete. 
m Pan ty $ "n vV39310. 2323.1. Mp MEE P V. E EAS ae) Pe ACRE AR tea 
J " H Ne 4 E NCMO CCP SUELOS eae 7, Uta DRACO TN, AV RR Pes ee ie CS RW ur wey ar. d 
dá CÁM E rd "7 re. 3 NH RE J A M2 Te oe ce 22 AL NAE ieee PAY, a r4 Aa COOPER E UTICA TODO Pea 
TP - F > i E E E E T E KU E TE T E E" A E X ME e AM REP Wn 2r ad ae wie AP TH gh 929-84 ROM X 211 LAUR AM OE LN nr De 
, p 4 i. 999 DE JC MN LEE TP P^ MAS nC it ne, c^ LATEIN Kr Pa STU RRA Seti PE 
LL: LECCE PR D NN C | E PR r RR a ed ri ONAN Pd n PE Y tum UE MS Y 
P H 4^ $$ Pour g^) y. 4. vef fte JURO nos dL iR VR fe CEU 8 ETE ia AGES. us n; curs 
UL L VERE 93 DULCE E: Ur EL £m w Duras Rasy Bee Paes DEDE Y D. ou An 
i H r" A X DURO UP osx 5 a ETO O 9 3 9, A raro i, Done: Ez M Ws e D LAS ahha fii Gy) Ant) 
- i L MU di SUA Kd X XD m meg. Xd A LAG FM E dy ) h^ * v SO REC dd * DE A " E PU. ^ a p FRE onda 
5 EN E f r a Th ae ee | EROS T2 ME k e wa " ^ UT: hi j^ D A. Y, DN cy A IY 3s B RL Ae ay ge f Et r1 
dé d > NN 8 LS a DE ) i H "MA Pu LC PETITS | neil i ex 4 He pt A. 4 iU RTA AS UD VL ig Tihs SCA KR Gag we 
E P i AERE S ee eet Ue ae ty Rs er oe Py Pas VONE ead tae date OF SERE Pr 
" $4? REPE ea ed ae DONEC MULCA Fed tee ee ECT ie UR iy Meet oer E "io ire 
; D "n Io’ 2a eT he Poe ee Pee k LSU ^fw RUE Me un Re ee oe EI A UL É 
" J oer, Fe EE WP $19 à 1 7 PUR EIE RE LRT: CAEN EVA oe En HOA SA UR Ari C on LT 
: QS x- a ^ 4 "H E AEA dry if ay yx) Va det RD cur vu Ponce Rv: y Fas wp, DITE A UE 
d SEL * uy f "fa W ni’ A EP "Yn PH E a DEP 2v KUETE Au pis CE Ya" "s Se iy Bt ge; d. hs 
. ny x A A TE? (EM ag +3 3 UE E FADA v^ D ID Wi 
ECT > Mo "y Pi ni PNE QC Py NS TEG) DON ee Ure rout A VY. DUE Pon Lets IE PUO 
D s EFA Si A XX DIC OUR US. LARA A UAE ETO ar es Do 43d Joi he VOR 
E u R d T * '* Ye ^ 1 dod - Pm A De . TL 
) Au^ v Y oua ds t ISIN ANS. "n. ur $ ^| ; £e dM US VE AV ded ae e EA 
b E sf $ IE * 4 Lun PUN US T | LA * t ». DR USE A. A f 
y 5 Tu UE " E K Wt kt yu y +> $. Ti VS fi D'Y Y, SU EA Y ARD IO d EA T unio y sie 
FX. ENCODED UU V OCC To JO A. yu UO ui PA ES T P DUE 
s r E Ww owe P db R TA + * yy ANLE Ez E MÀ beh 
AE PUES x n ZA Vi L aa k $33 X S ANSA E d V enn E DIN + ^h 
LE " “8 ; AT ) A Me A ee i ae mi bao fan n E ARIN d i Pet [Uo MES Pi AC Ear Y ges Aly V 
LI L] | a J * s " y A Ps F. A D 403 
, et AL PUE he J Me op 0 x BC EMEN UR d EE AAAA u A pU 
P E ? è à 3 à M, Y f NH KT me LE LUCR TE i 
: 7 de Me A " ^ vU H J? u y » (U^ OO Heg ^ tg vya Toda BOD PRY | CAA ee ee) Cre Pa A ce LV vb iu SEX deli tuy " 
EF'S ME OE PED TOLRE Ce a Tx een | Sus a DOM SA x os t [VV A QUAD aE PC 
r E Fy F . : DBAS, A wes ot bars Sota we Le TNR ¥ APS 4 WH TPF je " HA j^ Cea hee . NI Lu CUR. 
n LN IH er ie ae Ie ee ae PA Ji a LO "E OO mA MU UC OE eh ee eae ow yl vade g pd 9 $9 LOUPE ELE Aro 10.9 eo. Oy 
d E o UA ae OLR ACK ata ee ars OEC: oO era SS UIROS RUN OER Rie Cae agi gate EEG RC: 
P" j r K 4 S Oe oe M S MEN , v» 4 2. HE Ars KNEE HA 9.94 a v. LIE 6d H h 
E $ J agg Ce ILL? NITE qus M I JE RUE EN EC EL dors 
i hd 95529141 1,7 2) URL CE FAIDOSIDEC /y 4939; 4,9 4 Med bak DF nd 
JU S r CN Pa d UMP ar OEL9 ae CUT TC is Mee Oe ee eee ad te Ue hs J A Tod 
a 3 + $457 2.45 IA EA V) v.p 7175 0,1 CSL .0,9-% VE ata d 
EP 4 5 JUE TEE A 4 "M" f bed GR DK DOR Y 7 bet ee ae de a 
"ee t j MEL KL LAE PAR oh SOC 06 9o t^ 
ob y EEO ' PO 4 SS e uae. A t+ 9:0, 0 *7 Vin es £^ St . 
NS. ut s) F Te SEE ee Pe AE M UR PUE EET Ja E al i 
Sf 4.2) LE BRT 5 a 4 Wvfows. Me ODD Der D 
; r AES TEE PUN ZEE MR ? UR — p.v ghee D ye M Pea a a: NUM. X. a 
E 4 w A d 4v CEE DE F) E E AL ST Sae é «4 sf E TIPP RA TA wu SET 
P a h , Pathe pU CUR E QUE LE eur J4 FORTS N AULA 7 ENEA 
m2 v6 f " Hid TE EA AC * KA YUL M ar oet vy UN A Ye 
DNE. D MD a Pu B d dO M em LOEN E EN RAG pe UC ee A d i ft J Oy Lj Mt 
d ud LR e TAREE DAA ARR RAKAA OLUO kh 9M RN 6 m bep Ad ytd A IE 
i sà dd e E EA E E WM ES “tet UF te ee poe esr oe Li 
"^ Te S0 — PN d E ee Sethe ee a icy es 1X X SK na h FA AR) x E Lt ORE 
t pk p ^ a " "SF E TETE Eno "m. Vt "us "| T m, re wt. y ipid. 3 $;* eH LI le t 
5 * p LI n v . , iS » ? 
E 4 2 AL ^ JC NS me Or kee Nee I IY CRC ead be DEO OR UNA ; 
I EU be dat LAE EC 3x. r 528 3E ORDEN e UA ny OM . 
i FET E EEE AA em TAE ant CHREN A NER d s 
s s 5 P s LI MO N S o s AIVE E S ETA A AE OR R A t 4 betes 
[ELE T" é AEE 1 b ? Us C» 2x LL] JU IE LENT : v f à F 
i eee E pud yu i ! DAVIS RC y^ Vr Urb LULA ESO 
& 1 aa EX X A 7 oe rae ae Ff POLE G FE Palas O% 5,0 ie 
Fo .8 T EAE Y Us H De Y" URL da uA E 
^ ae ° ULIN Hus AE E M CP" PARA OR el Ake de ke 
P g UP : iy ONE CURL RELIANCE: ; 
UN g a AL ety J q oa pL ND RTI QUU 
MES OUT P. p de. à TARARE k EACT Y Nu P À UR 
ee ed 4 » y ED ' 4.8 ^s 4. LET j 9 4*9. t $9 EA AA Ni 
k J é a + 4b Re i Ft ty LA. LIAE Jo ub A i EL * 
a. i » H Q $ PA) i tis; at Oy E 2,453 M H NL A N y aey A 
V "t'.,5 $$ LE XN G SE IE r EO EAR d PE ar. LC EE. En e y $ Pd 
D , 4 ie te (ae ae © a LAN A D DR LP P PR AAA 
1 "Nan dor t oar, ^» 65i +f) vé 5 164 Lt P PARE D'S § 
BAR > I dL zs eb by ^ JUPE. BOE Vn - A AU n ie 
E e LI 3 PEE T EN í M EP E) LV. 3T Nun Ne P» i 
" adi T] LJ JUR SERRE ACED HERI t4 P o + I 4 
N 5 P Mr" " TP x A EI P Att V ^ D a! M JN 4% 
D E 0 E ¢ i XS MC F E at 
É f We. Ae hat a Die HE Le a 4 AP PERIE KAR, 
; ' AE JUL CA NR WA VW Me ^ sean et 
T4 Jd MIU Sang 14,7 a Ra Mii BuU e A Ur nq 
E A | 3 P M A 0 5 h d 
S n , PRT DOC E SCA C ELS ieee Leah taht OLN 
33 ‘ MEM L] Bue 1.3 Y, à 4 ry è 











NAVAL POSTGRADUATE SCHOOL 
Monterey, California 





THESIS 


A COMPARISON OF THREE NUMERICAL 
METHODS FOR 
UPDATING REGRESSIONS 


by 


Grigorios J. Raptis 


September, 1991 


Thesis Advisor: Dan C. Boger 
Thesis Co-advisor: William B. Gragg 





Approved for public release; distribution is unlimited. 





SECURITY CLASSIFICATION OF THIS PAGE. 


| REPORT DOCUMENTATION PAGE 


| 1a REPORT SECURITY CLASSIFICATION 1b RESTRICTIVE MARKINGS 
unclassified L.- 


2a SECURITY CLASSIFICATION AUTHORITY ! 3 DISTRIBUTION/AVAILABILITY OF REPORT 


2b DECLASSIFICATION/DOWNGRADING SCHEDULE 





Approved for public release; distribution is unlimited. 


4. PERFORMING ORGANIZATION REPORT NUMBER(S) 5 MONITORING ORGANIZATION REPORT NUMBER(S) 


6a. NAME OF PERFORMING ORGANIZATION 6b. OFFICE SYMBOL 7a NAME OF MONITORING ORGANIZATION 

Naval Postgraduate School (If applicable) 
OR 

6c ADDRESS (City, State, and ZIP Code) 7b ADDRESS (City, State, and ZIP Code) 

Monterey, CA 93943-5000 Monterey, CA 93943-5000 


Naval Postgraduate School 


8a NAME OF FUNDING/SPONSORING 8b OFFICE SYMBOL 9. PROCUREMENT INSTRUMENT IDENTIFICATION NUMBER 
ORGANIZATION (If applicable) 





8c ADDRESS (City, State, and ZIP Code) 10 SOURCE OF FUNDING NUMBERS 


= T NS Project No i ask NO m Work Umit Accession 
Number 
11 TITLE (Include Security Classification) 


A COMPARISON OF THREE NUMERICAL METHODS FOR UPDATING REGRESSIONS 


Program Element NO 





12. PERSONAL AUTHOR(S) Grigorios J. Raptis 














13a TYPE OF REPORT EFT ET TIME COVERED [ua DATE OF REPORT (year, month, day) E PAGE COUNT 
Master’s Thesis 1991, September 95 
16. SUPPLEMENTARY NOTATION 


The views expressed in this thesis are those of the author and do not reflect the official policy or position of the Department of Defense or the U.S. 
Government. 


17. COSATI CODES 18 SUBJECT TERMS (continue on reverse if necessary and identify by block number) 
ae GROUP SUBGROUP Hyperbolic Rotation, Gram-Schmidt Factorization, Reorthogonalization, Column 
oo Permutation, Updating 


19. ABSTRACT (continue on reverse if necessary and identify by block number) 


Three numerical procedures are presented for updating regressions. All three methods are based on QR factorization, but after that they use 
different philosophies to update the regression coefficients. Elden’s algorithm updates using only the triangular matrix R. This procedure does 
not use orthogonal transformations, but it uses hyperbolic rotations. The modified Gram-Schinidt QR process is used by Gragg-Leveque- 
Trangenstein’s method where the matrix with orthonormal! columns is stored and updated. Chan’s algorithm computes a column permutation TT 
and a QR factorization of a matrix A such that a rank deficiency of A will be revealed. Although the three methods are based on different ideas 
and can be used for different purposes their comparison shows that Chan’s algorithin is the only accurate one in the rank deficient case, and 

that Gragg-Leveque-Trangenstein's method is the cheapest and the most stable. 





20. DISTRIBUTION/AVAILABILITY OF ABSTRACT 21. ABSTRACT SECURITY CLASSIFICATION 
E] unciassirieo/uNuMimEo C] same as report C] omcusers Unclassified 
Dan C. Boger (408) 646-2607 OR/Bo 


DD FORM 1473, 84 MAR 83 APR edition may be used until exhausted SECURITY CLASSIFICATION OF THIS PAGE 
All other editions are obsolete Unclassified 


Approved for public release; distribution is unlimited. 


A Comparison of Three Numerical 
Methods for 
Updating Regressions 


by 
Grigorios J. Raptis 
Lieutenant Commander, Hellenic Navy 
B.S., Technical University of Athens 
B.S., Hellenic Naval Academy 


Submitted in partial fulfillment 
of the requirements for the degree of 


MASTER OF SCIENCE IN OPERATIONS RESEARCH 


from the 
NAVAL POSTGRADUATE SCHOOL 


Septembér 1991 i 


ABSTRACT 


Three numerical procedures are presented for updating regressions. All three 
methods are based on QR factorization, but after that they use different philosophies 
to update the regression coefficients. Elden’s algorithm updates using only the 
triangular matrix R. This procedure does not use orthogonal transformations, but it 
uses hyperbolic rotations. The modified Gram-Schmidt QR process is used by Gragg- 
Leveque-Trangenstein’s method where the matrix with orthonormal columns is stored 
and updated. Chan’s algorithm computes a column permutation II and a QR 
factorization of a matrix A such that a rank deficiency of A will be revealed. Although 
the three methods are based on different ideas and can be used for different purposes 
their comparison shows that Chan’s algorithm is the only accurate one in the rank 
deficient case, and that Gragg-Leveque-Trangenstein’s method is the cheapest and the 
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Ts INTRODUCTION 


A. BACKGROUND 

Regression analysis provides a variety of modeling 
techniques that allows the analyst to relate a dependent 
variable to one or more independent variables. The 
mathematical complexity of the model and the degree to which 
it is a realistic model will depend on how much is known about 
the process being studied and on the purpose of the modeling 
exercise. In preliminary studies of a process or in cases 
where prediction is the primary objective, the models will 
frequently fall in the class of models that are linear in the 
parameters. That is, the parameters enter the model as simple 
coefficients of the independent variables or functions of the 
independent variables. Such models will be referred to loosely 
as linear models. In many statistical problems, it is useful 
to express the dependent variable as a linear function of the 
independent variables. Furthermore, regression analysis can 
Summarize data and quantify the nature and strength of the 
relationships among variables. Regression analysis can also be 
used to predict new values of the dependent variables based on 


observed relationships. 


B. ANALYSIS OF REGRESSION 

As the reader will see in the following chapter, the main 
purpose of this thesis is to analyze the efficiency and 
numerical stability of some procedures for updating 
regressions. Most regression problems, however, require 
decisions on which variables to include in the model, the form 
the variables should take, for example, X, X^, 1/X, etc amd 
the functional form of the model. It is assumed that there is 
a set of t candidate variables, which presumably includes all 
relevant variables, from which a subset of r variables is to 
be chosen for the regression equation. The candidate variables 
may include different forms of the same basic variable, such 
as X and X^, and the selection process may include constraints 
on which variables are to be included. For example, X may be 
forced into the model if X* is in the selected subset. "Ehe 
a common constraint in building polynomial models. 

There are three distinct problem areas related to this 
general topic: 

B The theoretical effects of variable selection on the 
least squares regression results. 

B The computational methods for finding the "best" subset 
of variables for each subset size. 

B The choice of Sm ora (for the final model), or the 


"stopping rule". 


This thesis mainly discusses the second problem area of 
finding the more efficient and numerically stable 
computational methods. Also, we generally discuss the criteria 
for the "best" subset size choice, in other words the criteria 
for the "stopping rules". 

The simplest linear model involves only one independent 
variable and states that the true mean of the dependent 
variable changes at a constant rate as the value of the 
independent variable increases or decreases. Thus, the 
functional relationship between the true mean of Y,, E(Y;), and 


X; is the equation of a straight line, 


E(Y,) =B,+B, (X,) I 


where fj is the intercept, or the value of E(Y,;) when X - O, 
and f, is the slope of the line, or the rate of change in E(Y,) 
per unit change in X. The constants f, and f, are population 
parameters which are to be estimated from the sample of 
observations. 

The observations of the dependent variables, Y,, are 
assumed to be random observations from populations of random 
variables with the mean of each population given by E(Y,). The 
deviation of an observation Y, from its population mean E(Y,) 
is taken into account by adding a random error, €,, to give 


the statistical model 


Y,-p,*D,X,-*e,. (1.1) 
The subscript i indicates the particular observation unit, i 
— 1,2,...,n. The X, are the n observations of the independent 
variable and are assumed to be measured without error; that 
is, the observed values of X are assumed to be a set of known 
constants. The Y, and X, are paired observations; both are 
measured on every observational unit. The random errors €; 
are assumed to be normally, independently and identically 


distributed:[Ref. 1] 


€,-NID(0,0?). 


1. The Matrix Model 

Given the above regression model, we can write it in 
matrix form as Y = X6 + c, where Y is the n x 1 column vector 
of observations of the dependent variable. The n x p matrix X, 
where X — X;;7 ... = X, = 1, will be called the regression 
matrix, and the x;,;,'S are generally chosen so that the column 
of X are linearly independent; that is X has rank p. However, 
in some experimental design situations the elements of X are 
chosen to be O or 1, and the columns of X may be linearly 
dependent. In this case X is commonly called the design 
matrix. Also, B is the p x 1 vector of parameters to be 
estimated; and € is the n x 1 vector of random errors. 


Writing out the components of Y = XB + € yields 


Yi Xo M1 M2 +++ Hyp P. €i 














Y) X10 X231 X23 «5 Xp P. €z 

PE | N a (1.2) 
M la. X X2 m XD. Mi €En, 
(nx1) (ncc) (posl) (NT) 


The elements of a particular row, r, of X are the 
coefficients of the corresponding parameters in f which give 
E(Y,). The vectors Y and € are random vectors; the elements of 
these vectors are random variables. The matrix X is considered 
to be a matrix of known constants. The vector f is a vector of 
unknown constants and is to be estimated from the given data. 

Using the above assumptions on €,, the random vector 
€ has a multivariate normal distribution with a mean of the 
zero (0) vector. The variance-covariance matrix for any random 
vector of n elements is defined as an n x n symmetric matrix 
with the diagonal elements equal to the variance of €, and the 
off-diagonal elements equal to the covariance between €, and 


€ 





met Z be anx I vector of random variables Zj; Z2; Za; ==. 
,2,; the variance-covariance matrix of Z is the following 


mee n matrix. 


var(z,)  cov(z,,2,) ... cov(z,,z,) 


va cov(z;,,2,)  var(z,)  ... cov(z,,z,) (um 


cov(z,,2Z,) cov(z,,2,) ...  var(z,) 
The variance-covariance matrix of € is Io”, and since 
€ is independent and identically distributed, the distribution 


of € is written in shorthand notation as 


e~N(0, Io?), 


where I is the n x n identity matrix and c? is the common 
variance of all €,. Furthermore, since the elements of X and 
P are constants, the Xf term in the model is a set of 
constants being added to the vector of random errors, c. Thus, 
Y is a random vector with mean vector Xf and variance- 


covariance matrix Io: 
E(Y) =E(XP +e) =E( XB) +E(e) =xXB (1.4) 


Var (Y) =Var(XP+e) =Var (£) =Io? (1.5) 


Var(Y) is the same as Var(€) since adding a constant to a 
random variable does not change the variance. When € is 
multivariate normally distributed, Y is also multivariate 


normally distributed. Thus, 


Y~N(XB, Io?) (1.6) 


This result is based on the assumption that the linear model 


being used is the correct model [Ref. 2:.68]. 


2. The Normal Equations 

Before considering the problem of estimating 8, we 
note that we are referring to the model (2), where X;, is not 
necessarily constrained to be unity. In the case when X; # 1 
we have to use a notation in which i runs from O to p-1 rather 
than 1 to p. However, since the major application of the 
theory is to the case X;, = 1, it is convenient to separate f, 
Bron the other B.'s right from the outset. 

The oldest method of obtaining an estimate of f is the 
method of least squares. It dates back at least to Gauss and 
Exnsusts of minimizing »,€,¢ with respect to f (see Fig. 1); 
that is, setting gq = Xf, we minimize c'e! - |Y - qg|? subject 
to qcR[X] - W, where W is the range space of X (y:y-Xx for 
some x}. If we let q vary in W, |Y - qg|? (the square of the 
length of Y - q) will be a minimum for q - q' when (Y - q') is 
orthogonal with W (see Fig. 1). 


Thus 
x7(¥-q*) =0 aa) 


or 


1 We denote the transpose of a vector y by y!, and 


Similarly for a matrix transpose. 





Figure 1. The method of least squares consists of finding 
A such that AB is a minimum. 


Xo = xy (1.8) 
Here q is uniquely determined, being the unique orthogonal 


projection of Y onto W. Now given that the columns of X are 


linearly independent, there exists a unique vector (f) such 


that q' - Xf. Therefore substituting in (1.8) we have 


XX B =XY (1.9) 


the so-called normal equation(s). At this point we assume that 


X has rank p, so X!X is positive definite and therefore 


nonsingular. If so, equation (1.9) has a unique solution, 


namely, 


BP = (xX)" x'rY (1.10) 


Here ( is called the ordinary least squares estimate of f. 
Computational methods for calculating [ are given in the 
following chapters. Notice that [ can also be obtained by 
writing 
eTe=(Y-XB)7(Y-XB) = Y7Y-2B7X7Y+B7X xp 
using the fact that f'x'y - (f'x!v)! - y'’xf, and differentiating 
ete with respect to f. Thus from d¢€'€/0B8 = 0 we have 
-2X7Y+2X7XB = 0 
or 
X7XB = XY. 

This solution for f gives us a stationary value of €!e, and 
a simple algebraic identity confirms that [ is a minimum 
[Ref. 2]. 

The multiplication XTX generates a p' x p' matrix 

where the diagonal elements are the sums of squares of each of 
the independent variables and the off-diagonal elements are 


the sums of products between independent variables. The 


general form is 


n Y Xa ee 0 MX 
3$ Xa. eee, YO 5 => eee 
XX - Xa » X,:X75 Y Xi: Mire. NE i £i 


ly Xip 3o om We, PM y xi» 


Summation in all cases is over i - 1 to n, the n observations 
in the data. When only one independent variable is involved, 
X'X consists of only the upper-left 2 x 2 matrix. Inspection 


of the normal equations will reveal that the elements in this 


2 x 2 matrix are the coefficients of f, and fW. 


The elements of the matrix product X'Y are the sums of 
products between each independent variable in turn and the 


dependent variable: 


> Y; 
X XY, 


T = l2 
Moy yaa ( ) 
PE T, 
The first element, YY,, is the sum of products between 
vectors of ones (the first column of X) and Y. As is mentioned 


above, if only one independent variable is involved, XY 


consists of only the first two elements. 
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The fitted regression Y - xp, is denoted by 


¥Y (= {(Y,)]), and the elements of 
e = Y-Ŷ = y-xP = (21,-X(xX7X)? x")y = (1I,-P)Y, say, 


are called the residuals. The minimum value of ¢€'€, namely, 


ete = (Y-xB)7(y-xB) 


y Ty-2B7x TY+B7x xf 


y TY-B7x TY+P7 [x 7xB-x Ty] 


= yTy-B*xty = yTy-B7xTxB, (1.13) 


is called the residual sum of squares (RSS). Notice that Y and 


e are uniquely defined by f. 


C. VARIABLE SELECTION IN LEAST SQUARES ANALYSIS 

The purpose of the least squares analysis will influence 
the manner in which the model is constructed. Hocking [Ref. 4] 
relates six potential uses of regression equations given by 
Mallows [Ref. 3]: 

NH Providing a good description of the behavior of the 
response variable 

m Prediction of future responses and estimation of mean 


responses 
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@ Extrapolation, or prediction of responses outside the 
range of the data 

m Estimation of parameters 

m Control of a process by varying levels of input 

WB Developing realistic models of the process. 

Assume that the correct model involves t independent 
variables but that a subset of p variables, chosen randomly or 


on the basis of external information, is used in the 


regression equation. Let X, and f, denote submatrices of X and 


B that relate to the p selected variables. p. will denote the 


least squares estimates of f, obtained from the p-variate 
subset model and MS(Res,) the mean squares residual obtained 
from the p-variate subset model. After the above the following 
properties are summarized: 

W MS(Res,) is a positively biased estimate of 0^ unless the 
true regression coefficients for all deleted variables are 


Zero. 


a po is a biased estimate of f, and less variable Changer 





corresponding statistics obtained from the t-variate model 
[Ref. 4]. 

Stepwise regression methods are variable selection methods 
which identify good (although not necessarily the best) subset 
models, with considerably less computing than required for all 
possible regressions. The subset models are identified 


sequentially by adding or deleting, depending on the method, 
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the one variable that has the greatest impact on the residual 
sum of squares. These stepwise methods are not guaranteed to 
find the "best" subset for each subset size, and the results 


produced by different methods may not agree with each other. 


D. RESEARCH METHODOLOGY 
Given the regression model Y = X + €, where X is n xX p, 
a number of computational techniques have being suggested for 


solving the following steps: 


B Solve the normal equations Xx f - x'v. 


B Calculate the residual e = Y-xf. 


m Calculate the residual sum of squares RSS = e'e. 


m Update the regression model (that is, add or remove a 
row of X). 

B Add or remove a regressor (that is, add or remove a 
column of X.) 

B Calculate an  F-statistic for a general linear 
hypothesis. 

Solving the above steps is a common problem in a computer 
laboratory. These problems arise in a variety of areas and in 
a variety of contexts. Linear least squares problems are 
particularly difficult to solve because they frequently 


involve large quantities of data, and they are frequently ill- 
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conditioned*. The research methodology will be to describe 
acceptable procedures for updating regressions. The procedures 
to be chosen should be economically attractive and numerically 
accurate, in the case that the number of observations in the 
regression is large compared to the number of variables in the 
regression model. Each of these procedures (algorithms) has 


advantages and disadvantages which will be explained below. 


E. THESIS ORGANIZATION 
1. General 

The main purpose of this thesis is to bring to the 
attention of readers numerically acceptable procedures for 
adding and deleting rows and columns from a regression model. 
For the case of a single row to be inserted or deleted, the 
algorithms use simple techniques: plane or hyperbolic 
rotations and the Gram-Schmidt process. The algorithme for 
inserting rows are stable, but the problem of deleting a row 
may be ill-conditioned and some algorithms for this process 
may be numerically unstable. 

The need for updating regression results arises for 


various statistical or numerical reasons. When data are 


“A set of linear equations Bx = c is said to be ill- 
conditioned if small errors or variations in the elements of 
B and c have a large effect on the solution x. 
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arriving sequentially, it may be undesirable or impossible to 
wait for all the data before obtaining some regression 
results. In various time-series problems, one is interested in 
the changing relationships between variables. A regression 
model with a fixed number of lagged terms, as it moves over a 
series of data, generates a "window" on the sample, with a new 
observation added, and an old one deleted?, as the window 


moves to the next point in the series (Ref. 5]. 


2. Basic Schedule 

The procedures for updating regressions that we shall 
discuss here are the following: 

NH An algorithm based on the normal equations 
(Efroymson M. A., 1960[Ref. 6], Draper N. and Smith 
Ho 1966[Ref. 7]); this had been the standard 
introductory approach in regression courses. In this algorithm 
the regression coefficients are solved using Gauss-Jordan 
elimination on the normal equations. However, the problem of 
solving the normal equations is frequently much more ill- 
conditioned than the original problem of solving the 
overdetermined linear equations. 

m Stepwise regression analysis with orthogonal 


transformations (Lars Elden 1972) (REE «(8 |. In this 


?But in this case the problem can be ill-conditioned, 
because the rank can be decreased by this operation. 


L5 


algorithm a method is presented using QR-decomposition. In the 
QR-decomposition or factorization, we express a matrix as a 
product of an orthogonal matrix Q and an upper trapezoidal 
matrix R (actually R stands for right trapezoidal). Many of 
the ideas used in this method have been proposed by Golub 
(1965) [Ref. 9]. 

E The modified Gram-Schmidt triangular factorization 
(W. Gragg, R. Leveque, OV. Trangenstein 1978) 
(Ref. 10]. This approach is A = QR factorization with 
Q having orthonormal columns and R upper triangular. In this 
algorithm one is able to add or drop regressors (columns of A) 
or observations (rows of A) using essentially only the storage 
needed for A and to secure numerical accuracy possibly at the 
expense of additional computation and storage. 

E Rank revealing QR-factorization (T. Chan, Hansen 
1986). In this algorithm a method is presented for computing 
a column permutation, It, and a QR factorization, AI = QR, of 
an n by m (n > m) matrix A such that a possible rank 
deficiency? of A will be revealed in the triangular factor R 
having a small lower right block. Notice that a matrix A is 


rank deficient with rank deficiency d if it has at least d 


"I is a permutation matrix where [MeR™™ and is the product 
of P,P,...P,-;. Notice that P, denotes the matrix representation 
of the column interchange that precedes step i. 


"This means that the columns of the observation matrix A 
are not linearly independent. 


rG 


singular values. For matrices of low rank deficiency, the 
algorithm can reveal the rank of A, and the cost is only 
slightly more than the cost of one regular QR factorization. 
A posteriori upper and lower bounds on the singular values* 
of A can be used to infer the numerical rank of A. 

In Chapter II we discuss the general theory about 
updating regressions. A process for adding and dropping 
regressors follows. Techniques are presented with stepwise 
regression in mind, and we discuss how to compute the various 
quantities of statistical interest using the algorithms. 

Chapter III mainly covers computational details of the 
algorithms, which are used to compute the regression 
coefficients. In this chapter also a simulation is applied to 
the basic algorithms for validating the models by creating 
useful numerical results. In the same chapter we discuss the 
measures of effectiveness of each model based on the 
previously generated numerical results. Chapter IV concludes 


the thesis and offers recommendations. 


*PThe singular values of a matrix A are the "diagonal 
elements" of S, one of the three component matrices that A 
splitted to avoid singularity. 
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II. BACKGROUND THEORY 


A. GENERAL THEORY FOR FITTING A SPECIFIED REGRESSION 
Several techniques can be used to solve the problem of 
updating regressions. Efroymson's algorithm [Ref. 6] is one of 
the earlier ones used. It uses the method of Gauss-Jordan 
elimination on the normal equations A'Ax = A'b. A Mre 
efficient algorithm is to use the Cholesky factorization of 
the Gram matrix AČA into R'R where R is upper triangular. The 
solution to the original system is then found by a two step 


triangular solve process: 


Ry UTD Rx-y > A A-R Y AD. 


Another way to solve the above problem utilizes orthogonal 
transformations. This approach, called QR factorization, is 
the basis of the thesis algorithms and may be slightly slower 
than the normal equation approach, but is more stable 
numerically. The QR factorization of a matrix can be computed 
using the following methods: 

W Classical Gram-Schmidt (CGM). 

E Modified Gram-Schmidt (MGS). 

E Givens rotations. 


E Householder reflection (Elden's algorithm). 
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E QR with Column Pivoting for the Rank Deficiency case 
(Chan's algorithm). 

B Daniel-Gragg-Kaufman-Stewart method  (Gragg-Leveque- 
Trangenstein's algorithm). 
The latter algorithm uses rotations and the Gram-Schmidt 
process with reorthogonalization. 

Applying the QR technique on a square nonsingular system 
Ax = b we have (QR)x = b and by the associative law of matrix 
multiplication Q(Rx) = b. Defining y = Rx, it follows that Qy 
= b. Hence, Ax = b is solved in two steps: 

1. Solve Qy = b for the unknown y. 

2. Solve Rx = y for the unknown x. 
The second system is solved by back substitution. The first 
system is easy to solve since Q is orthogonal. Multiplying Qy 
= b by Q! yields QQ = Qb. Since Q is orthogonal, QQ = I and 
y = Qb. When A has more rows than columns and is of full rank 
the above process yields the solution of the least squares 
problem |b - Ax|,.7 = minimum. The key to the numerical 
stability is the orthonormality of the columns of Q. To avoid 
loss of orthogonality Gragg-Leveque-Trangenstein's algorithm 
applies reorthogonalization whenever |u|/|ul is small (say, 
Mess than. j/2/2), where u and ù is the vector to be 
orthogonalized and its orthogonal projection, respectively. 

Finally if A is rank deficient then the QR factorization 


need not give a basis for range(A). This problem can be 
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corrected, as mentioned below, by computing the QR 
factorization of the column pivoted version of A, AI = QR 
where II is a permutation, and applying Chan's rank revealing 
scheme. 

In the least square problem, statisticians and numerical 
analysts use different notations for the same entities, i.e., 
the matrix of the observations of independent variables, the 
vector of the dependent variables, the vector of random 
errors, etc. To avoid misunderstandings caused by using both 
notations in the following chapters, we built the following 
TABLE I of notational correspondents. 

Most numerical analysts use n to represent columns of a 
matrix and m to represent rows. Others interchanges m and n, 


and we will use this notation. 
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TABLE I 


NOTATIONS USED IN REGRESSION ANALYSIS 





E 


DESCRIPTION STATISTITIANS | N. ANALYSTS 


Independent var. matrix 


Dependent var. col. vector 


Vector of parameters to be 


estimated 


SS AB 
— 
Rows of indep. var. matrix com 5 
Cols of indep. var. matrix EL 
Vector of random errors mb 

redo 

NC 





Unique solution to the 


normal equation 


The vector of estimated 


means of the dep. var. 





aa: 
Paia! 
ee | 
pete] 
ETIN 


The idempotent n x n matrix 
The residuals vector 


Residual sum of squares 









SS (Res) 


Signific. level enter/stay pb 


The numerical analysts describe the vector p = Ax = b-r, 
as the "orthoprojector onto R(A)". 
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B. METHODS FOR UPDATING REGRESSIONS 

Before starting this section, let us introduce some 
terminology. The upper case letters denote matrices, lower 
case letters denote vectors, and Greek letters denote scalars. 


Thus we write 


€,) (Bs Pi) [m Quir Xu leu see 0, 
x-p4. ol, bzys els. |, AX a eee (2.1) 


\Ep/ \Bo) \Bp) \Np/ Goo Uu Mons cr eus) 


following the notational correspondents of TABLE I. For 


notational convenience we set x = B = b and the normal 
equations can now be written Xb = y or Ax = b. 


1. Lars Elden (1960) 

In Elden's algorithm the upper triangular matrix R in 
the QR-decomposition of A (in the model Ax = b) is determined 
by successive Householder transformations so that after k 
steps A has the form 

RD; qp Oo ae 
PLP... eee | 0 TE ai AUS 
(k) (1) (rye) 


R(? is upper triangular and T? is k x (m- k). Furthermore, 


T Jc. where W+ with [|w,,l, 1 is now chosen so 
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that k first rows are left unchanged and 


0 

E _ 440 

| A = Ae,, where |Al| = lanl. 
2) 221 


After m steps X will be completely triangularized: 


R 
Ep EAE 


In stepwise regression analysis, however, the final 


result may be of the form 


RP lp) 
OTA -| ) (2.2) 
0 A P» 
where RP is an upper triangular p x p matrix. The 


corresponding right hand side and the residual vector are: 


5 c oO O 
Q y = "I and r - Ds 


The regression coefficients are calculated from the 
system R(?b - c(?), In most cases the column vectors of X are 
not added to the subspace in the natural order, so that only 
after a permutation of columns is QA of the form shown in 


equation (2.2). 


a. Choice of Vector for Enlarging the Subspace. 
We shall add to the subspace the vector that will 


make the norm of the residual vector decrease as much as 
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possible. Since only the m - k last rovws of the matrix are 
changed in the kth step we can without loss of generality 
restrict ourselves to considering the first step of the 
analysis. We have the linear system of equations Ax - b, with 
A and b as in (2.1). Let a, 2 (04, 05, ..., 0,)' denob MEE 
jth column vector of A. We shall multiply (A: b) by an 
orthogonal matrix P = I - 2ww'!, where |w|;* — 1, so thats 5 
some j we have Pa; = @e,, where |A| = |a,l;. After this 
transformation the residual vector is given by the last n - 1 
components of Pb and since we were to decrease the norm of the 
residual vector as much as possible, j shall be chosen so that 
the first component (Pb), of Pb will be as large as possible. 
Now we have (Pb), —- (Pb)!e, - (Pb)!Ae,/A - b!P'Pa,/A - b!a,/A. 
Thus we shall select j so that the quantity 


(b7aj)* 


V, = 
J 
la;1,° 


will be as large as possible (Ref. 9]. 


b. Choice of Vector for Diminishing the Subspace. 


Suppose that after a number of steps (A: b) has the 


form 


*The Euclidean length of a vector is often denoted |xl;, 
also called the 2-norm, since the components of x are raised 


to the second power. 
|x|, = ER *XI4 -- -X. 
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(235 


0 A0 qth) 

where R“ is an upper triangular (K x k) matrix. Let r, = (Piy 
DNE 00090. 79) denote the pth column vector of Ry. 
Now if the jth column vector of A is removed we can by 
ENcceseijve rotations in the (j, j)*1), (j*1, j*2),...,(X-71,. K) 
planes transform R, to triangular form apart from the jth 
column. That is, we compute the coordinates of the remaining 
column vectors of A in the orthonormal basis where the jth 
base vector has been excluded. 


0 ees Ws 
Let Qj, denote the orthogonal rotation matrix in 


the (dj, j+1) plane which deviates from the unit matrix only in 
the elements qj; = Qj, jı = cos0, and qd,,4, - -dj4,; ^ sin8,. If 


0; gus chosen so that 


cos0,-— —Pi1 _ sin0,-——— Piu —— 


2 Z 2 2 
V Py, pea TP yea, p41 V P3, 1:1* P331, 11 
(6) (6) 
we get, Qj, 3115444 = Qj, jaa Ps, pur Py paar Pau pur Ore 70)* = 


/ 
= (pi, p.a se af Pyia, ys. Py,7+17979, s 9 8 , 0) sib where 


/ / 2 2 
Pj, j+1 = Py, 741 tP 441, 441° 


Further we have 


(04) / / 
Oy, jor Ij = (i, ,; sà (0 P34a,37 P3, 7) Pysi,g7 97 ss re) 


25 


where p!,, — cos0, pj; amd ptm — -sinG)9, ,. After k =m 


rotation we have transformed R, to the form 


Pia °» > © Pa, J- Pi,j Pi jaa ° * * Pale 
Pj-31,3-1  P3-1,3  Pyj-1,3:31 * * * Pyj-k 
/ / / 
Ro) demo eR Ora P3,J Pj,j«33 c v PUE 
a 0 a a 
/ 
; : '— Pk-i.k 
Px, 1 0 a 8 8 0 
(8-—] (0) 
where Q, = QUE fe 08 8 g¢ Qui. 


The preceding can be illustrated by the following schematic 


example (Kk = 5, j = 3) 


XXXXX XX 208 8 00 XXxXXxXxx 
Sie E A DN S) 

R = XXX = Xx = XXX 

XX te (lie x X 

x, | x + 1 
where x = any nonzero element, + = nonzero element being 

introduced, O = element being annihilated. 

To update the transposed inverse R7 = X! we 


multiply by the same rotation matrices from the left as in the 


following schematic example. 
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X x X 
x x rotation |x x WG) et tal OI See 
oM orn AM B) 
Ell Xx x m » x [] + - x X x 
xX xX Xx XXXX xcx Lux. 
xxrxxx XXX XX) ux x xxx 


Since Ra can be transformed to triangular form by permuting 


ix 


its columns, it is easy to see that (Rd will be triangular 


if the same permutation is performed on its rows. Therefore 
men th column of Q.(R^))* — (R) has only. one- nonzero 
EXEnent. Let 5 — (0;.—0, Epron €4) “denote the jth 
column vector of (R/?)T. Then Q,x,; = Ae. Since | x;l2 = | Q.x,| 
we have |A| = |x,ll2- The right hand side of (2.3) is multiplied 
by Q, and as the dimension of the subspace is decreased the 
increase of the norm of the residual vector will come from the 


Brae component (O,c*’), of Q,c"'. But 


T k Q, J Qx Q. J 


Thus in each step we have to find for which value of j 


v; (x; c (k)) 2 

Ix; 
is minimal, and if for this j the increase of the norm of the 
residual vector is not significant the jth column vector of X 
is removed from the subspace (the variable x, is deleted from 


the regression). 
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c. Inserting a Row (Observation). 
In some applications it may be desirable to insert 
a row to the matrix of observations after the analysis is 
finished and study what influence that row can have on the 
regression coefficients. This can be done without having to 
start from the beginning again using the QR-factorization that 
has already been done. The procedure can be illustrated by the 


following somewhat simplified example. Suppose we have the 


system of equations Rx = c, where 
Pia Piz © * > Pin Yi 
P22 * * * Pan Y2 
R= l f , c= 
P nn, Yn 


The system is augmented by a row (an observation): 


Pii Pio . Pin Yı 
Poo + © + Pan Y2 
Ro = r ; Co = 
P nn Yn 
Gr, €na;2 * * 5: WKCgauun pos 


If we multiply {Ry i5€9) bygga- rotation matrix Orar (0,) rns 


(1, n*1) plane where 8, is chosen so that 
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5 2 E / 2 2 
cos ei = Pai / pii*&yai,i , sin 8, DS pa,1 / Pii afii 


we obtain 


(1) (1) (1) 
pii Pio 5 5 — Pip 


P22 uon e Pon 


Ro” E Du R, = . . 
E 
0NES o MADE 
where 
pity = cos @, p, , + sin O, a,,, , j=1,...,n 
aha sein Op. > CED hi. J= 2an. 


After n successive rotations in the (1, n+1), (2, ntl), 


v e . , 


(n, nt1) planes 


(1) (1) (1) 
Dii Pi2 = = + Pip 


(2) (2) 
P22 >» + = Pan 


Ryo = One 82. Co RN (ODIERN- 


(n) 
Pon 


cy has been transformed: 


: a) TS) (n) p (n) 
cU? I (0) ,... go a (8.) ca = (Ci C2 rrey Cee Dens 
Dene have nehe n yeten R x =) (c,°%?,. 63°) 1o Ve eT 
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(the component b,4,U? belongs to the residual vector). 


d. Removing a Row (Observation). 

Sometimes it may be desirable to remove a row of 
the matrix of observations after the solution has been 
obtained. This can be done in a very simple manner without the 
analysis having to be done from the beginning again. Suppose 


that the matrix A of observation has been decomposed: 


R e 
T T 
A = ; b = (2.4) 
Qa | : Oa | A 
Let a, (vy = 1,...,m) denote the row vectors of X. Then the 
normal equations of the system Ax = b can be written in the 


following form: 


m m 
Y aʻa, = Mab. 


v=1 v=1 
If a, is the row of A that is going to be removed, let 


R 
s-[, | wheme j3UN— 4. 
P 


Then S! S =R R=a a R O IO RII IE 


Il 
f» 
< 
fy 
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Tt is easy to see that if the vector c in (2.4) is augmented 
by a component ib’, the corresponding right hand side of the 


normal equation will be 


m 
J a be 


S can be brought to triangular form by successive 


metr plications by matrices Timi: Tzar «~> 7 Inm Where T p is 
not unitary but complex orthogonal (T, mı = I). Consider the 
first step: 
Let 
Pii m 
1 0 
i oe (2:5) 
i 2 2 
VPii *pi 
0 T 
ia,, Pi 
Then 
/ / / 
Dii Pia += = = Pin 
P22 * * * Pan 
Ti, n+ S = s 
Pan 
/ / 
O fa,... . iapa 
where 
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1 _ (PaPa aptr) 


Piy z : doscdo ceret 
V P1i-7X*pi 
ila -p.,X, 4) 
a^, = p Pay Pii pi^ kg. o TS si 


2 2 
V P117 0p: 


After n transformations S has been brought to the form 
R! 
| P 
and the right hand side to the form 
c 
"I 


where we denote with double prime (") the final components of 


the normal equations. Then we have the system R"z = c" to 


solve. Notice that the component TD belongs to the residual 


vector. The hyperbolic rotations that are used to update the 
regression coefficient by working only on the matrix R are not 


unitary. Non-unitary transformations can destroy numerical 


Stability. 
2. Gragg-Leveque-Trangenstein (1978) 


These algorithms are based on the use of (orthogonal) 


plane rotations and the modified Gram-Schmidt process with 


32 


reorthogonalization. We compute the following quantities of 
interest in a regression analysis, with stepwise regression: 

B The least squares coefficients b, 

m The residuals r = y - Xb, 

mg The ANOVA table’, 

m The covariance matrix for the regression 
coefficients, and 

m Partial F-tests for the significance of a given 


variable in the regression. 


a. Stepwise Regression: Adding and Deleting 


Regressors 

Let us have entered k -.1 columns of X in the 
model y = Xb, and we wish to add the next column. This means 
that we want to find a k x k matrix Q, whose columns form an 
orthonormal basis for the span of the k columns of X, assuming 
that Q,, is known. The remaining columns of X are regressed 
upon the first k columns. Let R be a kx k triangular matrix 
which maps the orthonormal basis into the original basis. Also 
r, is the orthogonal projection of the vector y of dependent 
variable onto the space orthogonal to the space of the first 


k columns. In other words, r, is the residual vector. 


Mathematically, with k - 1 columns entered, we have factored: 


“Here, the residual sum of squares is computed directly, 
and not by subtracting the sum of squares due to regression 
from the total sum of squares. 


0 


Ry Tyas Cues 
X Y= Ox-1 Apr Tei] 0 I © (2.6) 
07 0 7 1 


where Q,., is n_x _(kK-1) and 
Qk-1 Qx- = I; (2.7) 
also the columns of A,, and r,, are orthogonal to the columns 


of Q,,. This means that 


Q,., 7A,.,70, Ons Ir, 470; (2.8) 


Also R,, is (X-1) x (k-1) upper triangular; and 


Ck-17Qk.i QUE (2.9) 


The above procedure was a step of the Gram-Schmidt 
process that can be viewed as a partial QR factorization. 
Although there are some similarities between the Gram-Schmidt 
factorization and Householder's factorization, there are also 
some important differences. We are going to discuss that in 
the next chapter during the comparison process. 

To update the above factorization (i.e., entering 
the kth column or, in other words, to add the kth variable), 
we perform the following steps: 


B Repartition: 
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Rees ty Cee, 
07 1 07 0 

Oo I 0 
o9 0 1 


Xoyv|- QNA, ri (2.10) 


m Normalization, and the equation (2.10) becomes: 


Rei Sk Ty, Ck, 


~ T T 
Qt. i a, A, Tta 0 led 0 O 
[24] 0 0 I 0 


070 07. ^T 


B Orthogonalization, with (dq, 


a,/ | a, | ) / to 
obtain: 
Rei Sk T, Ck-i 


r9 = oT a TA Tr 

Qi gx Aug (q.5A.) cu qur.) | | dx k dx k-1 
0 0 T 0 
0T 0 07 1 


B Relabel, so that the equation (2.10) can be 


written as: 


Ry Ty Cy 

O AFTO T 6 

007 1 
The above mathematical statements complete the 
insertion of the kth regressor. Furthermore the work required 
is essentially 2n(p - k) multiplications and additions. The 
storage space required and the part of storage affected by the 


algorithm is illustrated by these partitioned matrices. 


S 


Suppose now that we have to enter the kth variable 
(column) into the model. In this case the appropriate column 
to be added is the one which minimizes the norm of the 
residuals. We have that qi = 0 and from the 


orthogonalization and relabel steps, that 
Xe = yp My (Me Fp) OF Lyi = Fetlar Ty) Ay- 


By the Pythagorean theorem, 


|zilA 7 qm4^ 7 Cama (2 
so when we add the kth variable to the model the change in the 
residual is (qy'r,,)?. Recalling that qy s a,'/|la,|, the change 
can be expressed as (dqir,4)^ —- (a,'na/lal)^. So the chose 
column to add at the kth stage is a column a;, k<j<p, for 
which |a;'r,;|/|a,| is maximal. A test for significance should 
be satisfied for the above column to be entered into the 
regression. The work for choosing the column to add is on the 
order of 2n(p - k) multiplications and additions. 

To compute the regression coefficients we have to 
follow the next steps. A n x k matrix X, is formed by the 
first k columns of X. It has been factored X, —- Q,R,, where Q, 


is n x k with orthonormal columns, and R, is upper triangular. 





The projection of y onto the range of X can be written as: 


y-ry 7 Q,C, 7 Xb, - Q,R,b,. 
So the regression coefficient bh, can be computed easily by 


solving the following triangular system 
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RD, = Cy, 
requiring an order of 0.5k* multiplications and additions. 
The next step is the computation of the total sum 
of squares |yl?. A smart mathematical computation is taking 
place here. After computing the sum of squares due to the 


regression at the kth stage of the regression, as: 


Ebd? = OR byl? = RP? = jed? 

since c, is a k-vector, the residual sum of squares is 
computed directly as |r,|^, rather than by subtraction which 
would mean a loss of accuracy due to cancellation. This is 
recommended, especially if the mean residual sum of squares is 
to be used in sensitive tests for significance levels. The 
common criterion for terminating the selection process is the 
ratio of the reduction in residual sum of squares caused by 
the next candidate variable to be considered to the residual 
mean square from the model including that variable. This 
criterion can be expressed in terms of a critical "F-to-enter" 
or in terms of a critical "significance level to enter" (SLE), 
where F is the "F-test" of the partial sum of squares of the 
variable being considered. We can compute the "F-to-enter", 
Fe, for the kth variable by using the residual sum of squares 
before and after the insertion as in (2.6). Let RSS,., = 

(ry.1)'(ry,.1) and RSS, = (r,)'(r,) be respectively the residual 


sum of squares after each variable addition. Also let (n - m) 
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be the degrees of freedom for the last model, where n and m 
are the numbers of rows and the columns of the current matrix 
X of the independent variables. After the above, the "F-to- 
enter" is 


RSS,, - RSS, 


SE RSS,/ (n-m) 


To compute the covariance matrix of the regression 
coefficients b, using the normal equations, (X,'X,)7? usually is 
computed. For avoiding this numerically unstable process, a 
new procedure is followed in this algorithm. Using the above 


orthogonal decomposition, it is found that: 


(X R) aR D R Dy 


where R,«4 — (R/!)! - (Ru) l. The approach of the algermthneee 
to solve RV, = I, and compute (X,'X,)7 as V,'V,. Hence, V, is 


updated as follows: 


d li n 0 | ii : RTV 
T SPEEN V. ce 
2 OW Vie FM ag? Tag 


Thus, V, differs from V,, only in the bordering of a new row 
and column. The above updating of V, is mathematically 
equivalent to forwardsolving R,'V, = I directly. Notice that 


the computation of variances (i.e., the diagonal of the 
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covariance matrix) of the b, requires only on the order of k 
multiplications and additions. The variances of b, are needed 
for significance tests (i.e., t-test and F-test). 
Having entered k variables (columns) into the 
regression, we have the following factorization: 
Ry Ty, Cy 
Xy|-|Q.A,r,|o r O 
o1045 «1 
(KREIS 
where Q, is n x k with orthonormal columns, R, is k x k upper 


triangular, and 


QUI. zz Q,TA, zu, CK = OR. (2.12) 


Now, to delete the jth regressor, where 1 < j < k, the next 
steps are followed. The columns of X are permuted, so that for 


some permutation matrix P, 


Ry, (1,1) Re, (2,2) 843,331 Tk, Cka 


0 R 8 T c 

X viP- A r k, (2,2) 7J;:2 “kag “k2 
a nr 0 Oo I 0 (2.13) 

07 07 0 0 1 

(j-1) (k-j) (1) (p-k) (1) 


For k = 6 and j = 3, the right-hand factor looks like the 


following: 
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X X X XX X 
X X X X a 
Ey (1,1) Fr, ,3 Sga) xxx si 
O Rea 84,2) — xxx i 
XX 
: 


where x's denote possible nonzeros, and zeros are left blank. 
After that, by using a Givens transformation or (rotation), 
the subdiagonal of this matrix is annihilated. A Givens 


rotation is any orthogonal matrix of the form 


I 


G 
Jk = 
I rS Gea ke s ss s sj J+ 


1 
(3-1) (k-j-1) (n-k) (1) 


where G is of the form 


y O 0 o 
O 1 0 
0 i O 
Oo 0 " e a 0 —Yi 


The j and k subscript in G, correspond to the row numbers 
associated with the y's: The first y is in row k and the 
second y is in row j. Next the equation (2.13) is partitioned 


and after using the orthogonality of Givens rotations and 
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performing multiplications, the following equation is 


obtained: 


Ry (1,13. Fk, (33. 8531. Tei. Ck 
0 GR, (2,2) GS, 2 GT. 1 Gc, 2 
0 0 X 
07 07 0 07 1 


Oe, Osc Apa: (2:14) 


Furthermore (2.14) is repartitioned and relabeled to get its 


last shape: 
Rya Ta Cy-i 
Q,.1 Aka T-ıļ| 0 L 0 
07 07 1 


(k-1) (p-k+1) (1) (k-1) (p-kK+1) (1) 


Note that (2.14) holds with k replaced by k - 1 and that the 
deleted column is ready to reenter the regression at any time. 
The above algorithm uses approximately n(p+2k-33) 
multiplications and additions. 

However, the chosen column to drop is that one 
which yields the smallest increase in the residuals. This 
increase for the jth independent variable is |8,|? |v,l?, where 
B, is the jth regression coefficient, ana |v; is the jth 
diagonal element of (X,'X,) !|. Thus the column to be dropped, if 
a test for lack of significance is satisfied, is the column j 


for which |8,|/|v,| is smallest. 
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b. Adding and Deleting Observations 
Now n rows of X have entered into the regression, 
and an additional row x’ and observation n of y is going to be 
inserted. In this situation, a factorization has taken place 


and we have: 


Xa+ Yn| = = xT n 
xT n oF 1 O}\oT 1 


where Q, is n x p with orthonormal columns, R, is p x p upper 
triangular, also 

OnO = Deand pier Creen = 0 mea CeO vat (2.15) 
Furthermore the Givens rotations are applied and after 


relabeling, repartitioning and performing some elementary 


operations we get the following: 


07. gi 





Raw Cnet 
Qu ftYq ' 


Notice that (2.15) holds with n replaced by n -* 1 and that the 
above algorithm requires n + 1 additional storage locations 


for q. However, we can easily update V - R', since 
Vs 
I-R,TV,-(R,T x) PX 


and we apply the Givens rotations to get: 
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T 
n+l Vi i 





(RIT | 


V 
I = (R,7 x) 6% j 
0 


Now we want to delete the jth observation. Then if 
Ime the jth row of X, Nn is the jth entry of y, q' is the jth 


row of Q, and p is the jth entry of r, a permutation matrix P 


can be obtained such that: 








X Vn-1 R, Cy 
P Xx, TES -|P|Q, Ta 
xe wm 0e 
KC 
Ó $ R, C, 5 R n n 
- T oT 1 a Q 2l o" 
q Pp q7 1 po^ 1 


Next, we apply the Gram-Schmidt process. After this we have: 


1 
&) =0 = ( 1 =q Ta : / oG=-Oq, Y-qd'£, Ta-i -f-Yd; 


when o = O, (@7,@) is chosen by a special feature of the 


orthogonalization code. Now if we use this orthogonalization, 


we get: 
I q O|R, Ca R, €, 
2a 5:9 4d F107 o Yio" 6 /| =| 2 @ z,,||07 $*|. (2.16) 
x" qn q7To jp AO'O 1/07 1 a'o po i1 


Next we choose Givens transformations such that, (q', o)@ = 
(q7, W)Gp pi- - -Gip = (07, 7). Since [ (d', o)| = 1, it follows 


that 7 = +1. So the equation (2.16) can be written as: 
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GT 0ol/G OVR, C, 


um Xnaa| 7 9 q Ia o2 Y 
x TEE q7 o p AO" 1A0" 1/407 1 
Ro. Cpga 
= |Q; 0 F,4|l sT y |- (252020) 
0T « p 07 1 


because the matrix 





T “lar = Qi d 
q^ o 0T « 


has orthonormal columns it follows q' = 0 and because we have: 
R, E R, B Roi 
G 07 = Gy pe s Gp, pia 07 xm - s 
S 


where R,, 1S an upper triangular matrix. However from the 


equation (2.17) above, we obtain the desired factorization 


which is: 


RUP e 
Xa nA Qua Zn-1)\ 07 1 


and we can recover the dropped row by: 


xT = t87, N= Ty + p. 


Finally, we have to update V = R™. Now, 
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Since 7 = +1, it follows that uv" = 0. Hence R,,'V,, - I, and 


V,., is the desired update [Ref. 10]. 


3. Tony F. Chan (1986) 

An algorithm is presented for computing a column 
permutation |l and a QR factorization AI = QR of an n bym (n 
> m) matrix A such that a possible rank deficiency of A will 
be revealed in the triangular factor R having a small lower 
right block. For matrices of low rank deficiency, the 
algorithm reveals the rank of A, and the cost is only slightly 
more than the cost of one regular QR factorization. An upper 
and lower bound on the singular values of A are stated. These 
can be used to infer the numerical rank of A. 

A very useful factorization of an n by m (n > m) 
matrix is the QR factorization, given by AI = QR, where [eR™ 
is a permutation matrix, QeR™™ has orthogonal columns and 
satisfies Q'O - I,, and RER™ is upper triangular. 

If A has full rank, then R is nonsingular. In many 
applications in which A is nearly rank deficient, it is 


desirable to select the permutation so that the rank 
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deficiency is exhibited in R having a small lower right block. 


For if R is partioned as 
zs 2 
0 R4 
where R; is d by d, then it is easy to show that 
O -z1 (<4) < [&4];, where we have used the notation 0, to denote 
the ith singular value of A, with 0, 2? 0, È ... 2? O 
Therefore, if JR,,], is small, then A has at least d small 
singular values and thus is close to being rank d deficient. 


The converse, unfortunately, is not true. In other words, if 


A has d small singular values, then it is not guaranteed that 


a given QR factorization of A has a small |R,l. Leta,e 32» 


be the matrix of order n illustrated below: 


T ce © -C 

0 1 EC B s zC 

. E a . . . —C 
A -odiaco(l sic EM) ; 

O0 . 0 0 NI 


where s and c satisfy sf + cf = 1. For n= 50, c= 0.2, we have 
o,(A,) = 10%. On the other hand, A, is its own QR factorization 


and obviously has no small R, block for any value of d. 
Besides being able to reveal rank deficiency of A, a 
QR factorization with a small R,, block is very useful in many 


applications, such as in the rank deficient least squares 
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problem and in the subset selection problem. Therefore a 
variety of techniques have been proposed to compute it. Since 
OR factorization is essentially unique once the permutation II 
is fixed, these techniques all amount to finding an 
appropriate column permutation of A. Perhaps the best known of 
these is the column pivoting strategy. Although this strategy 
is usually very effective in producing a triangular factor R 
with small |R.||, very little is known in theory about its 
behavior, and it can fail on some matrices. Chan's algorithm 
does not require computing the SVD of A and is most closely 
related CO one recently developed by Foster 


meer. 11). 


a. Revealing Rank One Deficiency 


Assume that A is nearly rank one deficient. We 
would like to find a column permutation of A such that the 
resulting QR factorization has a small pivotal element r,,. It 
turns out that this permutation can be found by inspecting the 
Size of the elements of the singular vector of A corresponding 
to the smallest singular value 0,. This procedure was first 
pointed out in 1976 [Ref. 12]. 

Assume that there is a vector x e HR? with |x|, = 1 
such that lAx|, = €, and let I be a permutation such that if 


Ix = y, then |7,| = lyl. and Iyl = |xll, = 1. Such a x can be 
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obtained by the LINPACK? condition estimator. Now, because 


of these we have |n| > /1/n and furthermore 


Q'Ax - Q'AIL Ix - Ry -| ` 


Poon 


Therefore, 
e = [axl = IQ "Axh = IR yh > Paal. 
from which we have the result, Ip,,] s /n £. Now let v e P^ with [v,] = 1 


be the right singular vector of A corresponding to the 


smallest singular value 0,. Then we have 


[Av], =o 


A 
Therefore, by the above, if we define the permutation Il by 
|v) |a = vi. 

then All has a QR factorization with a pivot p,, at least as 
small as no, in absolute value. 

Since only the permutation IJ is needed, it is not 
necessary to compute the SVD of A in order to find v exactly. 
In practice, one can use a few steps of inverse iteration to 


compute an approximation to v from which the permutation II can 


be determined. In the more interesting case where 0,«9,,, the 


Together, LINPACK and EISPACK represent the state of the 
art in software for matrix computation. 
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inverse iteration should converge rapidly. This suggests a two 
pass algorithm in which one first computes any QR 
factorization of A, then performs inverse iteration with R to 
find an approximate v, then determines II, and then computes 


mae OR factorization of All. 


b. Revealing Higher Dimensional Rank Deficiency. 


In this section, we consider the case where A is 
nearly d deficient, with d > 1. Our goal is to find a 


permutation [| such that if 


is the QR factorization of AI, with R;; € R"?, then |£,l is 


small in some norm. 

A natural way to extend the one dimensional result 
of subsection a is to repeatedly apply the one dimensional 
ameeerithm, for d = 32,2,...,Ra4, the leading principal 
triangular part of R. Suppose that we have already isolated a 
small d x d block R,,. To isolate a small (d+1) x (d+1) block, 


we compute, using the one dimensional algorithm given in 


subsection a, a permutation P such that R,P = QR, is the QR 


factorization of R,P and where the (n-d, n-d)th element of R,, 


is small. Then with 
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n=n(? z 
0 I 


where, P is a permutation matrix such that P e R, 


|\(P7v),| = IP7vl, and the singular vector v € R! corresponding 


to Onin(Ri1) with |ull, = 1, and 6, 2» O,&,(Ri)), also 
Dm 
0 Ij 


After the above, it can be easily verified that 


Afi 


5 P Q. , 
0. RE 


is the OR factorization of All. 


To make the above procedure more understandable, 
the updating process is illustrated for n = 5 and i = 2. The 
permutation is defined by moving the second column of R to the 


last column. Thus 


NC XN XX 
X X X X 
gH,-| lx x *, 
OD x + 
E 
where again, x = any nonzero element, + = nonzero element 


being introduced and O = element being annihilated. So working 


from the left in planes (1,1+1), (i+1,1+2), ... , {(n-1,nj)eye 
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annihilate the subdiagonal of RII, by rotations. Simultaneously 
we compute Q - QQ,. The computation f = IP, is trivial, like 


that of RIl,. If 0,:,(Ri,) 1s "tiny" then we have for n= 5, 


X X Xx X X 

AL X X 

S X X X 
A= OR fl’ =O fi? 

K X 


Qc 


with € also tiny. If € is negligible we may drop the last row 


EN and the last column of QO to get 


x xxW LL. 
Xe Le - 

xxū 

x O, 


with © having orthonormal columns. Now if we annihilate the 


last column of the upper trapezoidal matrix, as before, from 


the bottom to top, and then drop the last columns of the [I and 
R matrices, we obtain A = OR fl". 

The above algorithm, to produce the desired QR 
factorization, is based on the following two assertions for i 
EN PEUn-1, ucesch-dtl: 

8 R, has a small singular value so that the (i,i)th 


element of R,, is guaranteed to be small. 


o 


" The last [i.e., the (n-i)th] row of Q/R, is 


small. 

If these two assertions are true, then the lower 
(n-i-1) x (n-i-1) block of R is small and we have the desired 
QR factorization. But these two assertions are true because of 
the following lemmas proved in (Ref. 13]. 

Let B € R™ be a matrix containing any subset of k 


columns of A. Then 


O.1,(B) = 0,(B) < o,(A). 


Also if the matrix W = [w,-.41,---,W,] € R™* has been 
computed by the above algorithm then it should satisfy the 


following properties: 


e |w = 1, 
= (w); = O for j >i; 
= | (wi)il = Iwi la>=1//i, 


= JjAllw,|]; — 6, € 0;(A). 

Finally, Chan's algorithm [Ref. 13] 
computes a permutation IJ and a QR factorization of A given by 
All = QR where the elements of the lower d x d upper triangular 


block of R satisfies 


ji 
|r| « 96,73 * y 203-9 o,/k « 207? o,/n 
k=1 
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fore T-G < 1 < j < i. 

Very often, it is desirable to be able to infer the 
rank of A from its QR factorization by estimating the small 
Singular values of A from the triangular factor R 
(Ref. 14]. For this, we state bounds on the singular 
values of A in terms of the matrices R and W, given by the 


following inequality, for 1 <j < d, 


O a-4-+1 


-1 
ee O a a S Rah e 9; 710082) h. 
JJIKO 7 | 
2 2 


where, W € R™, R € R™ and computed by the rank revealing 
algorithm, R;? and W; denote the lower right j by j upper 
triangular blocks of R and W respectively. These bounds are 
proved in [Ref. 13]. 

Now using the obtained All = QR factorization we can solve 
the least square problem Ax - b as follows. Let r be the rank 
of the n x m matrix A, then the regression coefficients are 
given by: 

xXx = Ilo (RO 1) 30545) 
mere 1 = 1, 2, 3, ... , r.[Ref. 15] 

The work for the All = QR factorization is given by: 

W(r) = m*(n - m/3) + Imr + 2m’r, where I is the number of 
inverse iterations used at each step. Usually I = 2 is 


sufficient in practice. So 
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W(r) - m^(n - m/3) + 4m’r. (2.18) 


C. PROCESS 

In the problem of finding a vector x c R" such that Ax - 
b where the data matrix A c R"" and the observation vector 
b € R” are given and m > n when there are more equations than 
unknowns, we say that the system Ax = b is overdetermined. 
Usually an overdetermined system has no exact solution since 
b must be an element of range (A), a proper subspace of R”. 

This suggests that we strive to minimize |Ax — bI EENSM 
Suitable choice of p. Different norms render different optimum 
solutions. However, much progress has been made in this area, 
and there are several good techniques available for 1-norm and 
o-norm minimization. The oldest method for solving the full 
rank least squares problem is the method of normal equations. 
The accuracy of the computed normal equations solution depends 
on the square of the condition number, and the algorithm is 
not always accurate. 

For the above reason some techniques are established 
based on the QR factorization method. The Householder and 
Gram-Schmidt QR approaches to the least squares problem are 
more stable than the normal equation method. Furthermore 
techniques for rank revealing QR-factorization (i.e., Chan's 
algorithm discussed in subsection 3) can give a solution to 
the subset selection problem even if A is nearly rank 


deficient. 
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In the following chapters we will try to compare the 
results of the above algorithms computationally to verify the 


effectiveness of each on the least squares problem. 
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III. RESULTS AND COMPARISONS 


A. PREVIOUS COMPARISONS 

There have been a number of studies comparing Householder 
and Modified Gram-Schmidt QR techniques to form an orthonormal 
basis of R(A) [Ref. 16]. Moreover, in  Elden's paper 
there is a comparison between his and Efroymson's algorithm 
for updating regressions. There are comparisons in work and 
accuracy between the Classical and Modified Gram-Schmidt 
methods in the Gragg-Leveque-Trangenstein paper. Tony Chan 
also compares his column permutation QR factorization 
algorithm with the regular QR algorithm, in case of required 
work  [Ref. 17]. Chan's comparison is applied only on 
the QR factorization part of the regression updating 
procedure. 

Furthermore, an algorithm for solution of the subset 
selection problem is presented in a technical support package 
of NASA Tech Briefs [Ref. 18] that can be mentioned 
as a modified Chan's algorithm. So a comparison can be done 
between these algorithms. Finally a justification of the use 
of the reorthogonalization in Gram-Schmidt QR factorization is 
presented in the Daniel-Gragg-Kaufman-Stewart paper 


[Ref. 19]. 
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B. MODEL AND DATA 

The data in TABLE II-1 will be used to illustrate the 
model selection methods in the full rank case. In the table 
below there are two matrices, the 13 x 5 matrix X consisting 
of a column of ones, followed by the 4 column vectors of the 
observations on the independent variables, and the 13 x 1 
column vector Y of observations of the dependent variable. In 
order to conform with previous results, the data utilized in 
this experiment is identical to that used in Elden's algorithm 
[Ref. 8]. We began the stepwise regression analysis by 
inserting and deleting columns of the observation matrix X and 
we continued by inserting and deleting rows of the model. 

The data in TABLE II-2 and TABLE II-3, have one and two 
dependent columns respectively. We obtained the first of them 
by subtracting columns 2 and 3 of the data in TABLE II-1 and 


the second by repeating the first column of the same data. 
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| | TABLE II-1 


CASE STUDY 1 (INPUT DATA 


MATRIX OF OBSERVATIONS 


COLUMN 
VECTOR OF 
MATRIX OF OBSERV. ON THE INDEPENDENT VARIABLES | OBSER. ON 
DEPENDENT 
VARIABLES 


26.00 | 6.00 | 4i.00 | 78.50 
29.00 | 15.00 | 52.00 | 74.30 


.00 | 52.00 —- 
| 1.00 | 7.00 | 52.00 | 6.00 | 33.00 | 95.90 - 
| i.o0 | 11.00 | 55.00 | 9.00 | 22.00 | 109.20 - 
| 1.00 | 3.00 | 71.00 | 17.00 | 6.00 | 102.70 - 
| 1.00 | 21.00 | i.00 | 4.00 | 26.00 | 115.90 _ 
| 1.00  ! 1.00 | 40.00 | 29.00 | 40.00 | 3.80 
11.00 | 66.00 | 9.00 | 12.00 | 113.30 - 


1.00 t "To. 00 68.00 8.00 T220 109.40 





To obtain the data in the Table II-2, (a non-full rank 
matrix), we replaced the column 4 of the Table II-1 with a new 
column, column 6. This new column was obtained by subtracting 
columns 2 and 3 of the matrix of observations in Table II-1. 
Now the new matrix of observations has a linearly dependent 


column and so it is rank one deficient. 
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TABLE II-2 


CASE STUDY 2a (INPUT DATA) 


MATRIX OF OBSERVATIONS 


COLUMN 
VECTOR OF 
MATRIX OF OBSERV. ON THE INDEPENDENT VARIABLES | OBSER. ON 


DEPENDENT 
VARIABLES 


| 1.00 | 6.00 | 3.00 | 71.00 | -68.00 | 102.70 | 
| 1.00 | éo.o0 | 1.00 | 31.00 | -30.00 | 73.50 — 
1.00 115.90 

| 1.00 | 34.00 | 1.00 | 47.00 | -39.00 | 83.90 — 

66.00 


00 12700 10.00 68.00 =o OO 109.40 





Furthermore in Table II-3 we inserted once more the column 
1 of the matrix in Table II-2 as a new column 7. This column 
took the third place in the matrix and reduced its rank to 


rank two deficient. 
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TABLE U-3 


CASE STUDY 2b (INPUT DATA) 
MATRIX OF OBSERVATIONS 








COLUMN 
VECTORZOF 
MATRIX OF OBSERV. ON THE INDEPENDENT VARIABLES | OBSER. ON 
DEPENDENT 
VARIABLES 


| 1.00 | 6.00| i.00| 3.00 | 71.00| -68.00| 102.70 - 


26.00 | 1.00| 21.00 -26.00| 115.90 


33.00 52.00| -45.00| 95.90 
34.00 40.00| -39.00| 93.10 


12.00 | 1.00. 11.00 66.00| -55.00 113.30 





C. SIMULATION 

As mentioned, the three compared algorithms use a QR 
factorization technique to solve the regression problem. But 
the three techniques used have different philosophies and 
different advantages and disadvantages. However, to obtain the 


numerical results needed for the comparison we used three 
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analogous sets of MATLAB codes. It is important that Elden's 
algorithm is not a "full" updating regression algorithm 
because there is no way to update  regressions using 
Householder QR factorization. Elden adds or deletes columns or 
rows at the beginning or the end of the observations matrix 
and does not explicitly update the  Householder QR 
factorization. Because of this we did not implement Elden's 
algorithm. Instead, we used the basic QR factorization 
algorithm to solve each individual problem. For column 
insertions and deletions on the right of the matrix our 
algorithms are computationally equivalent with  Elden's 


algorithms. 


D. VALIDATION 

A numerical example was needed to illustrate the theory 
developed in Chapter II, so for each of the three algorithms 
we ran the corresponding MATLAB codes on the data of the case 
study. All computations are done on a 286 PC in double 
precision, with a relative machine precision of about 107°. 
Notice that during the stepwise procedure on case study 1, 
there was no rank deficient case, so there was no illustration 
of the effectiveness of Chan's algorithm to reveal the 


numerical rank of a given matrix. 


61 


E. MEASURES OF EFFECTIVENESS 

1. Case Study 1 (Full Rank Case) 

In order to compare the performances of the algorithms, 
the STSC'S Statistics/Graphics package STATGRAPHICS 4.0 was 
used on the same data. The results of the algorithms and 
STATGRAPHICS stepwise variable selection procedure are 
presented in the following tables. 

STEPWISE REGRESSION ANALYSIS WITH ORTHOGONAL 
TRANSFORMATIONS 


NUMBER OF OBSERVATIONS 13 


STEP N- 1 


TABLE III-1. Summary Statistics for Stepwise Selection of 
Variables for the Data in TABLE II after step 1. 


VARIABLE ENTERED AFTER THE COLUMN OF ONES: Column 5 


(2) (3) 
GRAGG- 
ELDEN LEVEQUE- CHAN 
TRANGENST. 


| CONST. | 147.5679312 | 117.56/9212341177.56/9312 | TIAS ieee 


VARIABLE 
COEFFICIENTS 


ALPHAS į -0.7381618 =O. 7 3o tole -0073816 E3 =0. 738 T65 


RES. 

SUM 883.8669169 | 883.8669169 | 333.3669169 | 881.89616 
SQUAR. 

Fe 22.1985202 22. 798353202 22.79805202 22.7985 
LEVEL 





2018 22.7995203 2 2/99/52 0) 22 . 7995203 22:80 
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STEP N 2 
TABLE III-2. Summary Statistics for Stepwise Selection of 
Variables for the Data in TABLE II after step 2. 


VARIABLE ENTERED: Column 2 


Mp nue 2) Le c6) ee 


CONST. | 103.0973816 103.0973816 Į 120340973816 103.097382 


VARIABLE 
COEFFICIENTS 


ALPHA5 20 6139520 =0 269639536 -0.6139536 *0.615954 
ALPHA2 154399593 1.4399583 1:399593 1.439958 


RSS 74.76219722 34. 762358122 74.7621122 74.7621 
Fe | 108.2239093 | 108.2239093 | 108.2239093 108.22399 


E GEN. | 176.6269531 1 4656269253189 M 1 46-.62695J.1 176.6270 
| MODEL 





STEP N 3 
TABLE III-3. Summary Statistics for Stepwise Selection of 
Variables for the Data in TABLE II after step 3. 


VARIABLE ENTERED : Column 3 


(4) 
CONST. || 71.6483069 | 71.6483069 | 71.6483069 | 71.648307 


| VARIABLE 


i 
i 


COEFFICIENTS 


ALPHA5 -0.2365402 720102565 27 -O8 2369402 -0.236540 
PHA2 T. €5T9379 1.459379 11.4519379 1.451938 


ALPHA3 0.4161098 0.4161098 0.4161098 0.41611 


47.9727294 47.9727294 47.9727294 47.9727 
5.0258647 5.0258647 5.028647 5.0259 


166.8316801 1/166.8516801 7 166.8316801 166.8317 





63 


STEP N 4 
TABLE III-4. Summary Statistics for Stepwise Selection of 
Variables for the Data in TABLE II after step 4. 


VARIABLE REMOVED: Column 5 


EE EE eee 


CONST 52.57753489 5235775239 52.5773489 52 «D> le Jee 


VARIABLE 
COEFFICIENTS 


ALPHA2 1.4683057 1.4683057 1.4683057 1.468306 
ALPHA3 O.36622505 0 -6622505 0.6622505 0.662250 
| 


RSS 57.9044832 57.9044832 57.9044832 57.9085 | 
Fe | 1.86328 73 1.8632 873 1 : 8632873 1.8632 


-F GEN. || 229.5036971 | 229.5036971 | 229.5036971 | 229.5040 
| MODEL 





STEP NM 5 
TABLE III-5. Summary Statistics for Stepwise Selection of 
Variables for the Data in TABLE II after step 5. 


OBSERVATION ENTERED: Row 3, (a second time) 


(2) 


- E m n — -— ——— mm p — oto 


: rs - P SL z E 
CONSIL.XI] 52.6817201 Do 46 09 7 2404. 924'600 "209 52.680201 


VARIABLE 
COEFFICIENTS 


ALPHA3 
Rss | 59.9550974| 59.9550974 | 59.9550974 | 59.9551 —- 
re | -—— | -—— [| —— | -—— | 


I. ESGENSIII 250.34377701]8/250.3437770 M1250,5437770 |N1250 7528 
MODEL 
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STEP N 6 
TABLE III-6. Summary Statistics for Stepwise Selection of 
Variables for the Data in TABLE II after step 6. 


OBSERVATION ENTERED: Row 2, (a second time) 
CONST. 53.039011? 53.03 80932 5303801 Y2 53.038011 


VARIABLE 
COEFFICIENTS 


ALPHA2 1.4484905 1.4484905 1.4484905 1.448490 
eee d 0.6549147 0.6549147 0.6549147 0.654915 
| 60.8055442 60.8055442 60.8055442 60 “6055 





SmeP N" 7 
TABLE III-7. Summary Statistics for Stepwise Selection of 
Variables for the Data in TABLE II after step 7. 


OBSERVATION REMOVED: Row 1 
ret ee ee a 
| CONST. | 53.0380116 53.0380116 3 3.0 380716 53.038012 


VARIABLE 
COEFFICIENTS 


ALPHA2 1.4484905 1.4484905 1.4484905 1.448490 
ALPHA3 0.6549147 0.6549147 0.6549147 0.654915 


EE presens 60.8055442 | 60.8055442 | 60.8055442 | 60.8055 


912.7948771 | 322.7948772) į 312.7948771 | 342.7949 
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H 












2. Case Study 2 (Rank Deficient Case) 

Here we want to test the thesis algorithms on a rank 
deficient case. So we used the singular matrices in Tables II- 
2 and II-3. The first matrix has one linearly dependent column 
(column 5) and the second two linearly dependent columns 
(column 3 and 6). We applied the corresponding algorithms' 
MATLAB codes to those matrices. The results of the algorithms' 
effectiveness are shown in the following tables. 


TABLE IV-1. Summary statistics for the rank deficient case, 
(case study 2a). 


RESULTS OF CASE STUDY 2a REGRESSION 


(1) (2) (3) 
GRAGG- 
ELDEN LEVEQUE- CHAN 
TRANGENST. 


CONST. | 3.2155*10!?? | -146.000000 | 71.6483069 


VARIABLE 
COLFETCTENTS 


-0.0353*10!* 2.000000 20.2365402 


| -0.0194*10!* 2.5676*10!^ 1.4519379 
-0.0173*10!*? -2.5676*10!^ 0.4161098 
-0.0132*10!? | -2.5676x*10!" 0.0000000 


5.7494*10°° 45.8653934 47.9727294 


9.2803*107?* | 116.4231889 ! 166.8316801 
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TABLE IV-2. Summary statistics for the rank deficient case, 
(case study 2b). 


RESULTS OF CASE STUDY 2b REGRESSION 
"m. (1) (4) 
CONST. 5.786*10"" -7.6870*10" 71.6483097 


VARIABLE 
COEFFICIENTS 


ALPHA5 || -0.007242 3.8663*10!^| -0.2365402 
ALPHA7 || -5.786*10** 7.6870*10*> 0.0000000 


ALPHA2 0.708*10!* -3.8663*10!^ 1.4519379 


ALPHA3 | -0.708*10!* 3.8663*10!^ 0.4161098 
ALPHA6 || -0.708*10!" 3.8663*10!^ 0.0000000 


59.17750 a7 Ser Ons 47.9727294 


212799352033 3.953*107°! 166.8316801 


A discussion of the above results is given in the 





following chapter. 
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IV. CONCLUSIONS AND RECOMMENDATIONS 


A. GENERAL 

The major emphasis of this thesis was to examine the 
performance of three algorithms for updating regressions. As 
mentioned, the algorithms examined were Elden (Householder), 
Gragg-Leveque-Trangenstein and Chan. After running the 
corresponding MATLAB codes we obtain the results shown in 
Tables III and IV. The discussion of these results will be 
divided into two categories: the accuracy and stability, and 


the number of computations. 


1. Accuracy and Stability 
a. Full Rank Case 
The results of the case study 1, in Table III, show 
us that no one algorithm uniformly outperforms any other with 
regard to accuracy in the full rank case. All algorithms give 


exactly the same results in all steps. 


b. Deficient Rank Case 
Comparing the results of the case study 2 (TABLES 
IV-1 and IV-2) with the results of the case study 1 step 3 
(TABLE III-3) we can see that they are identical except for 


the coefficients of variables ALPHA6 and ALPHA7. These 
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variables are used only in the case study 2 and they 
correspond to the linearly dependent columns of the matrix of 
observations (see TABLES II-1 and II-2). The values of these 
variables' coefficients in both rank deficient cases are zero. 
Because of the above we can say that the results of Chan's 
algorithm in the rank deficient case without the linearly 
dependent columns are the same as the other algorithms' 
results in the full rank case. However in case study 2 (rank 
one and rank two deficient cases) only Chan's algorithm gives 
reasonable results. The results of the other two algorithms 
are totally wrong. Now it is easy to understand that the 
philosophy of Chan's algorithm is to drop the linearly 
dependent columns of the matrix of observations and work with 


the rest of them to obtain the regression coefficients. 


2. The Number of Computations 
TO compare the volume of computation of the present 
algorithms we consider that all m column vectors of ann x m 


matrix A are added to the subspace. 


a. Computing Regression Coefficients at Each Step 
In the case where the regression coefficients are 


computed at each step the work of algorithms is as follows. 
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The number of flops’! of Elden's method is about 


2 


2nm^, (2(n - m)m? + 4m°/3 for producing the OR factorization, 





m?/3 for updating the inverse of the triangular matrix and m?/3 
for computing the regression coefficients). Elden's algorithm, 
as mentioned in Chapter II subsection Bl, is updating 
regression without using at each step the matrix Q but working 
only on the triangular matrix R. This procedure avoids 
computations on matrix Q and so is cheap. 

If the data matrix is not very rank degenerate Chan 
uses same work as Elden on the updating regression procedure, 
(Householder without updating Q). Using the equation (2.18) 
for this case with r - m we have a total work for the rank 
revealing QR factorization of W(m) = nm’ + 11m?/3. Also the 
amount of flops for computing the regression coefficients is 
m?/3. In other words the work for Chan's QR factorization of 
RII is about nm? + 4m? flops. 

Gragg-Leveque-Trangenstein's algorithm requires 8nm 
flops for column insertion, O flops for column deletion, 3m(m 
+ 2n) for row insertion, m(3m + 14n) flops for row deletion. 
In this case the regression coefficients are computed at each 
step we have to add the work of solving the triangular system 
Rub, = c,. That requires a total of m°/2 flops. This amount of 


flops is because, as mentioned in Chapter II subsection B2, 


Flops number means the whole number of multiplications 
and additions. 
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this algorithm updates the regression using the matrix Q for 
keeping the numerical stability. That means that in the worst 
case the algorithm requires a total work of about 6m’ + 28nm 
+m?/2 flops. 

After the above we can see that Elden's and Gragg- 
Leveque-Trangenstein's algorithms outperforms Chan's 
algorithm. As mentioned, in this thesis we always examine the 
case in which the observation matrix has more rows than 
columns (i.e., n > m). Keeping a ratio n/m = 2 Gragg-Leveque- 
Trangenstein's algorithm is the cheapest for n bigger than 18. 
This means that for big applications this algorithm is the 


cheapest. 


b. Computing Regression Coefficients Once 


In the case where the regression coefficients are 
computed once we have a lower order work for this computation 
which can be ignored. So the total work for the algorithms is: 
Elden's 2nm^ - m°/3, Gragg-Leveque-Trangenstein's 6m* + 24nm 
and Chan's nm? + 11m°/3. Keeping the same ratio as above Gragg- 
Leveque-Trangenstein's algorithm is still the cheapest for n 


bigger than 15. 


B. RECOMMENDATIONS 
In this thesis we examined three procedures of QR 


factorization to update the variable selection problem. The 


pA 


above factorization can be also used on several other 
mathematical and statistical procedures. Important research 
can be done on using the QR updating algorithms to solve the 
linear programming problem using the simplex method and 
Karmarkar's algorithm. 

As mentioned in Chapter III section D, we used Householder 
codes to simulate the procedure of Elden's algorithm. A closer 
comparison can be done by a full simulation of this algorithm. 
Finally, it should be interesting to extend this thesis in 
case where we have an observation matrix with fewer rows than 


columns (izes en eama 


Ho 


APPENDIX A. MATLAB CODES FOR HOUSEHOLDER'S (ELDEN) ALGORITHM 


Dies codes Sqrit.m,, Grhup.m, dqrucd.m, qruri.m, dqrurd.m 


qrorth.m and rot.m were reproduced here with Professor's Gragg 


permission. 


o? o9 og? o? o? o? o9 oe 


o? o9 


oo 


oe 


oo 


oe 


function W = qrhf(A) 
W = qrhf(A): 


W is a pack formed array which contains the HOUSEHOLDER QR 
FACTORIZATION of A. We have A = QR with Q unitary and R 
upper trapezoidal with nonnegative diagonal elements. The 
commands R = W, (Q R] » qrhup(R), executed by qrh, unpack. 
However this is not necessary, and in fact it's inefficient, 
for most applications. 


Copyright (c) 16 February 1991 by Bill Gragg. All rights 
reserved. 


qrhf calls sgn. 


begin qrhf 
{n m] = size(A); 
for k = 1- min(m, n) 
q = k:n; u = A(q,K); Pinot); t = u(l); 
lf 9:220 
t = u(1); s = - sgn(t); u(1) = t - r*s; 
t = abs(t); t sol tatr, d = sqrt(2*t); 
A(q,K) = u/d; da = r*žsqrt(t); u = u/d; 
Ik m 
P- kti: nm: W — A(q,p); W-W - u*(u'*W); 
pc WC S) Ve = Se, A(q,p) = W; 
end 
end 
end 
W = A; 
end qrhf 


function [Q,R] -» qrhup(R) 


[Q R) = qrhup(R): 
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o9 o9 o9 o9 o9 o9 o9 o o9 o9 o9 oo 


o 


oe 


oe 


o9 o9 oe? o9 o9 o0 


oo o? oo 


oo 


Computes the ELEMENTWISE QR factorization of A given the 
output R. Thus [Q R]=qrh(A) :=qrhup(qrhf(A)) is essentially 
equivalent with matlab's function qr, except that we execute 
an inexpensive unitary diagonal scaling to make the diagonal 
elements of R nonnegative. It is NOT USEFUL to have Q in 
elementwise form to solve the LS problem, norm(b - Ax) = 
minimum. The purpose of matlab having Q in such form may be 
to avoid having to explain how the Householder least squares 
algorithm really works. 


Copyright 16 February 1991 by Bill Gragg. All rights 
reserved. 


qrhup calls sgn. 


begin qrhup 
(n mJ] = size(R); S = - sgn(diag(R)); e = ones(n-m,1); 
Q = diag([s;e]); q = m-*1:n; z = zeros(n-m, 1); 
EOOt2 - sqrt(2); 
for k = min(m,n):-1:1 


q= [kaq]; u = R(q,k); r = norm(u); 
ift-me 0 
u = u/(r/root2); T = Q(q,dq); Q(dq,dq) = Tuam QN E 
end 
R(K,k) » r; 
Pk aon 
R(kK+1:n,kK) = 2; 
end 
ZzZ - [2;0]; 
end 
end qrhup 
function W = sgn(Z) 


w = sgn(Z) or W - sgn(Z2): 


For z a complex number w is 2Z/abs(z) if Z = 0 and + 1 if 
z = 0. Thus sgn is the same as matlab's sign function EXCEPT 
when z = 0. We always have abs(w) = 1. W is the elementwise 


(Schur) sgn function of the complex matrix Z. 


Copyright (c) 19 January 1991 by Bill Gragg. All rights 
reserved. 


sgn calls no extrinsic functions. 


begin sgn 
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oo 


of of oO 


o9 o9 o9 o9 o9 o9 o9 o9 o9 o9 oe 


oe 


oo 


oo 


oo 


O° oO 


W = sign(Z); p = find(Z == 0); n = length (p); 
W(p) = ones(n,1); 
end sgn 


Total flops (scalar case, see csgn): 
Real case: O flops. Complex case: 1 sqrt + 4 mults + 1 
add. 


function [x,r] = qrhlsl(W,b) 
fe 2) = gqrhisi(W,b): 


Given that the HOUSEHOLDER QR FACTORIZATION of A, is stored 
in packed form in the array W, qhrs SOLVES the least squares 
problem norm(b - Ax) = minimum for x. It also efficiently 
computes the residual r := (b - Ax). We assume that rank(A) 
= m <= n, where A has m columns and n rows. 


qhrs calls no extrinsic functions. 


begin qrhs 
[n m] » size(W); root2 - sqrt(2); 

Forward solving: computing Q'b = H(m)H(m-1)...H(1)b. 
for k = 1:m 


q = k:n; u = W(q,k); rl = norm(u); d(k) = ri; 
gf m2 
ú = u/(rl/root2); v = b(q)?) b(q) = y —^u*(u*'*v); 
end 
end 
Backsolving: solving Rx = b1 := b(1:m) for xX. 
s = - sgn(diag(W)); x = zeros(m,1); 


x(m) = s(m)'*b(m)/d(m); 
for k = m-1:-1:1 
PKD], X(K) s) Sb Sj Wile pie x(p)y7 aCe), 


end 

Computing the residual r. 
if m <n r= b —- W*x; else r = 0; 
end 

end qrhls1i 


function [Eb,RSS,F,Fe] - regres2(Q,R,P,RSSp,b,r) 


[Eb RSS F Fe]2» regres2(Q,R,P,RSSp,b,r): 
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Eb = E(b), 1s the expectation of beta of the linear least 
Squares problem, norm(b - Ax) = minimum, RSS is the residual 
sum of squares and F the statistic for a general linear 
hypothesis. Also Fe is the significance level to enter, SLE, 
"F-to-enter", for the new variable, on where A - QRP' is the 
QR general FACTORIZATION of A. p = A*x is the orthogonal 
projection of b onto R(A) and the residual vector r - b - pr 
—b - Ax, is the orthogonal projection of b onto N(A'), the 
orthogonal complement of R(A). It is assumed that A is not 
a zero matrix. 


regres2 calls rows, cols, ones and norm. 


o9 o? o9 o? o o9 o9 o o9 o? oQ o? o9 o? 


begin regres2 


Eb = r * Q*R*P'*x; 

A = O*R*P'; n = rows(A); m = cols(A); 

e — ones(n,1); RSS - r'*r; v = Q'*b; 

f = (norm(v)-abs(e'*b)/sqrt(n))* 
(norm(v)+abs(e'*b)/sqrt(n)); 

df = (n-m)/(m-1); F = (df*f)/RSS; 

Fe = (abs(RSSp-RSS)*(n-m))/RSS; 


9 


% end regres2 
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APPENDIX B. MATLAB CODES FOR GRAGG-LEVEQUE-TRANGENSTEIN'S 
ALGORITHM 


function ([(Q,R] = qruci(Q,R,a,i) 
kR] = qgruci(0,R,a,1): 
UPDATES the QR factorization A = ỌR when a is INSERTED as 


COLUMN i of A. Q has orthonormal columns and R is upper 
trapezoidal with at least as many columns as rows. 


o9 o9 o9? o9 oe 


Copyright (c) 20 July 1991 by Bill Gragg. All rights 
reserved, 


oP oe 


oe 


qruci calls qrorth and rot. 


oO 


begin qruci 


(r m] = size(R); n = length(a); (q s t] = 
darorth(Q,a)? 
R(:,i+1:m+1) = R(:,1:m); R(:,0)9 9 87 m=mi+ l: 
MI rmn 

Q = [Q q]; R = (R;zeros(1,m)]; 

r= r+ 1; Ria) = t: 
end 
for k = r-1:-1:i 

P= kil:m; d -= etal, force 


rot (R(q,1)):; 
R(q,1) = [(t;0]; 
Q(:,qd) *G; 
R(€a, P) = G'*R(q, P); 

end 

IC 39997 
q = itl:r; ad = diag(R); d = 
sgn(d(q)); 
D = diag(d); OT = Oad) =); R(q,:) 
D'*R(q,:); 

end 

$ end qruci 


Q 
Il 


[C - s';s C'];  Q(:,q) 


fünction [9,R; a] = qrucda(Q; Ri) 


% [Q R a] = qrucd(Q,R,1): 


d. 


oo o9 oo o9 o9 o9 oe 


oo 


oo 


oo 


OP oO A oO oO 


oe o9 o? oe 


UPDATES the QR factorization A - QR when COLUMN i is DELETED 
from A. Q has orthonormal columns and R is upper trapezoidal 
with at least as many columns as rows. The deleted column is 
called a [Ref 19]. 

Copyright (c) 20 July 1991 by Bill Gragg. All rights 
reserved, 


qgruca calls Tot. 


begin qrucd 


if nargout > 2 a = Q*R(:,i); end, 
[r m] = size(R); R(= a = ie m =m- 1; 
for k = i:r-1 
p =*k+1:m; Cp eed [c s t] = rot(R(q,K)); 
R(q,k) = [t70]; G = [C -95';S CIE 
QU: nd). m Qi AG 
Pek << qu R(q,p) = G'*R(q,p); end 
end 
it m 
GS eb q = 1:r; 
O = QUU R = R(q,:); 
else 
S = sgn(R(r,r)); Q(*,r) — Qs SE 
R(r;:) "= S'IR(r e): 
end 


end qrucd 


function Dove Orid 9R, 279) 
[Q R] = qruri(Q,R,a,3): 


UPDATES the factorization A = QR when a' is INSERTED as ROW 
j of A. Q has orthonormal columns and R is upper trapezoidal 
with at least as many columns as rows. 


Copyright (c) 20 July 1991 by Bill Gragg. All rights 
reserved, 

qruri calls rot. 

begin qruri 


m = length(a); {n r] = size(Q); 
ea a; O(c n;:) -= [zeros (l1,r):; QO(j:n-1, E 
pou: T Q(:,r) = zeros(n,1); Q(j, r) = 1; 
R - [R;a']; 
for Kk psp] 

p = kł+1:m; = [k rJ]; 

[c S t] = rot(R(q,k)); 

Rita, ik) = [eo G = [c -s';s c!']; 

o D OE G 
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oo 


o? o9 o? o? o9 o? 


o? o9 o9 


oo 


oe 


oo 


iek <m R(d P) =~ C'ZR(d; p): end 


end 
fe. ro >-m 
r=re-l; d Soir. 
DNO SH R — R(q,:); 
else 
s = sgn(R(r,r)); Q(:, r) = Q(:,r)*ž*s; 
Rui) SUXRO S); 
end 


end qruri 


function ([(Q,R,a] = qrurd(Q,R,j) 
[OQ R a] » qrurd(Q,R,j): 


UPDATES the QR factorization A = QR when ROW j is DELETED 
from A. Q has orthonormal columns and R is upper trapezoidal 
with at least as many columns as rows. a' = A(j,:) = Q(j,:)R 
is the deleted row. 


gaurd calls grorth and rot. 
Copyright (c) 20 July 1991 by Bill Gragg. All rights 
reserved, 


begin qrurd 


[n r] » size(Q); [r m] » size(R); Eure] Qa etm 
Ie rsen 
q = zeros(n,1); GC coms r=r + 1l; 
OC: r) =- grorth(0,qd):? R(r,:) = zeros(1,m); 
end 
for k » r-1:-1:1 
p= kim ; q = [kr]; 
[c S u] = rot(Q(j,q)); 
Q(J),r) = u; = [-s c';c s']; 
Oera = Ce, 4) 4G, 
R(q,p) = G'*R(q,p); 
end 
r=ro- 1; q = 1:r; ODO Eo ys: 
a = R(rt1;,:)';. R= R(q,;:?:): 
da = diag(R); d = sgn(d); D = diag(d); Q = Q*D; 
R = D'*R; 


end qrurd 


function lga r,S,k] qrorth(Q,a) 


(qr s k] = qrorth(Q,a) or q = grorth(Q): 
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o9 oo o9 o9 o9 o o9 o9 o? o? o9 oo o9 o9 o9 oe 


oe o9 oe 


oo 


È 


© 


ORTHONORMALIZES column a against the columns of Q using 
REORTHOGONALIZATION. It is assumed that Q has (nearly) 


orthonormal columns. Let [n m] = size(Q). An error message 
occurs if m» n. For m <= m we have 
[Qa] - [OQ g)[I r] 
[ s] with Q'q= 0, s >= 0 


and, if m< n; 
q'q = 1. 
If m = n we take q := 0 and s := 0. If a is not input we 


take a := 0. For m < n and a = O we get a unit vector q in 
the orthogonal complement of the range of Q. 


qrorth calls no extrinsic functions. 
Copyright (c) 20 July 1991 by Bill Gragg. All rights 
reserved, 


begin qrorth 


[n m] 7 size(Q); 
if nargin < 2 a = Zeros(n,1); end 
if m» n error('Q has too many columns.'), end 
if m == n q = zeros(n,1); r= Q'*a; s = 0; k = 0; 
return, end 
norma = norm(a); t = norma/2; 
r = Q'*a; b = a - Q*žr; s = norm(b) ; KIRI, 
if norma == zflag = 1; k=0; end 
while 0 == 0 
if s > t q = b/s; break, end 
if k > 4 error('Process did not terminate in 5 
iterations.'), end 
if s <= t*eps/10 
[u j] = min(norms(Q')); 
if s == 
s = eps/4; b(jJ) = s; 
else 
Exp b(3) + s+eps/4; 
end 
norma = Ss; k= Kk + 1/2; 
end 
t = 5/2; b = b — Q*(Q'*b); s = norm(b); k=k+ 1; 
end 
end 
if zflag r = zeros(m,1l); s = 0; end 
end qrorth 


80 


o0 o9 o9? o9 o9 o9 o9 o9 oO o? o9 o9 o9 o9 o9 oo o9 oo 


oo oe oe 


oo 


oo 


oo o9 oo 


iíunction [c,9S.r]  rot(x;,y) 


(C S rj 


TOE(X, Y) Or [C S Y] 


rort(z 


J: 


Carefully computes the Gram-Schmidt QR factorization 


Note that 


Li Ge a 


[x] 
[Y] 


[x] 


[Y] 
[C e 
[s c' 


[c] r 


[S] 


] [x] =: QR 


] [9] 


is a full OR factorization. This is a "tool of the trade" in 


computational 


Copyright (c) 28 October 1990 by 


reserved. 


linear 
ROTATIONS and that 


[ c' 


[7S 


algebra. 


Note that 0Q 


s'] 
c] 


bot calls no extrinsic functions. 


begin rot 
if nargin < 2 


C = sign(x); 
ifyo 
if X> y 
t = y/x; 
u = 1/u; 
else 
t - X/y; 
u = t/u; 
end 
else 
r = x; if 
end 
end rot 


Total flops: 
Real case: 
Complex case: 


Y 
S 


l 
al 


X2) IX 
sign(y); 


sqrt(1 + 
u*c; 


sqrt(l + 
u*c; 


sqrt + 


SCC 
x 


CXC); 


txt); 


end 


[x] 
[Y] 


Bill Gragg. 


[r] 
[0] 


; end 


abs (x); y 


x*u; 
v*s; 


E 


I y*u; 
S v*s; 


5 mults + 1 add. 


sqrt + 11 muakts + 3 adds. 
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and  Q' 


are 


All rights 


= abs(y); 
v = t/u; 
v = 1l/u; 


oo 


o? o? o? o9 oo 


oo o? o9 oo 


oo 


oo 


function [x, r] = qrmgsls(O RD) 
[x r] = grmgsls(O R bD): 


Solves the norm(b - Ax) = minimum for x given the QR 
factorization A = QR and computes the minimal norm rnorm 
using the MODIFIED GRAM-SCHMIDT process (mgs). 


qrmgss calls gfsb. 


begin qrmgsls 

Compute the "Fourier coefficients" c = Q'b and the residual 
vector r= b - Ax = b - QORX = bd - Oc = (I - QQ')b, WITHOUT 
COMPUTING x. 


r = b; 
for i = 1:cols(Q) 
q= ool) t = q'*r; c= (ce 
end 
Compute r and backsolve for x. 
x = gfsb(R,c):? r= b - Qc; 


end qrmgsls 
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Q9 o9 o9 o9 o9 o9 o9 o9 o9 o9 oo oO o9 o9 o9 o9 o? o9 o9 o9 o9 o9 o9 o9 o9 o9 o? o? 


oo 


oP oP 


oo ge 


APPENDIX C. MATLAB CODES FOR CHAN'S ALGORITHM 


function (Q,R,Pi,rank,W,delta) = rrgqr(A,tol) 
R = rrqr(A,tol) 
(QO,R,Pi,rank,W,delta] = rrqr(A,tol) 


Compute a rank revealing QR factorization of A: 
PoP OsRe— Oper 1) R12 J], 
[ 0 R 22 ] 
such that Q is m-by-n, R is triangular n-by-n, and 

1) R 11 is rank-by-rank and well-conditioned, 

2) rank - numerical rank of A with respect to tol, 
defined as the number of singular values greater than 
tol, 

3) norm(R 22) is of the same order of magnitude as the 
(rank+1)'th singular value. 


Also, return a matrix W whose columns are orthonormal and 
span an aproximation to the null "N" space of A, and the 6 
delta(rank:n1:2) containing lower and upper bounds for the 
last n-rank-*1 singular values of A (the first rank-1 rows of 
delta are zero). 


If no tol is specified, sqrt(n)*norm(A,1)*eps is default. 


This program is an implementation of the algorithm described 


in the paper: T. Chan, "Rank Revealin R factorizations", 
Lin. Alg. Appl. 88/89 (1987), 67-82. 


The use of the factor c max ratio was suggested in: C. H. 


Bischof, & P. C. Hansen, "Structure preserving and rank 
revealing QR-factorizations", SISSC, to appear 

(m,n] 2» size(A); delta - zeros(n,2); c max ratio - 10; 

if (nargin--1), tol = sqrt(n)*norm(A,1)*eps; end 


if (nargout>4), W = [({]; end 


Compute an initial QR factorization A*Pi = Q*R. 
OR, Pi) = GriA)> O = O(:,1en)? R = R(lin,:); 


Prepare for the iterations. Estimate smallest singular 
value Of R. 

nu = 0; i= n; 

(sest,v] - ccvl(R); if (nargout»5), delta(n,1) = sest; end 
Loop until a singular value estimate larger than tol is 
found. 

while sest « tol, nu = nutl; 


Update the matrix W. 
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oo o9 o9 


oe? o9 o? oe 


oo 


oo oe 


oe 


o9 o9 o9 o9 oo oe 


if (nargout>4), W = [[v;zeros(nu-1,1)],W]; end 


Find the element in v with greatest index, numerically 
within a factor c max ratio of the numerically largest 
element. 
vmax - norm(v,inf); 
for k=1:-1:1 

if (vmax «- abs(v(k))*c max ratio), break, end 
end 


If necessary, generate the permutation that brings the pivot 
element of v to the last position, apply the permutation to 
W, Pi, and R, compute a new QR factorization of R(1:i,1:i), 
and update Q and R. 
if (k-i) 
p= [kt1:i;k]; PI UKk T) 
if (nargout- 4 BU 1) 
for j=k:i-1 
[es] -ogensgi(E() EE m 
[R(J,J:n) ,R(3J*1,3:n)] - 
app g left(c,s,R(j,j:n) ,R(3*1,j:n)) ; 
R(j+1,j) = 0; 
[Q(:,j) Q(:,)+1)] = app_g_right(c,s,Q(:,)) O(I TEE 
end 
end 


Pi(:,p)? R(:,K:1) = R(:7PIE 
W(p,:); end 


Provide an upper bound for the i'th singular value. 
if (nargout>5), dđelta(i,2) = norm(R(i:n,i:n)); end 


Estimate the smallest singular value of R(1:1-1,1:1-1), 
which is a lower bound for the (1-1) 'th singular value. 
i = i-l; [sest,v] = ccvl(R(1:i,1:i)); 

if (nargout>5), delta(i,1) = sest; end 


end 

Finish the computation. If nargout < 2, return R. 

if (nargout>5), delta(i,2) = norm(R(rank:n,rank:n)); end 
rank = n - nu; if (nargout>4), W = Pi*W; end 


if (nargout < 2), Q = R; end 
This algorithm is described by T. Chan & P. Hansen [Ref 15]. 


function [x,r] = basic(Q,R,Pi,b,rank) 
[x r] » basic(Q,R,Pi,b,rank) 
Compute the basic solution. If the RRQR of A is 


A*Pi — O*R —00 [ RSIE Reale ie, 
( 0 R22] 
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o? o? oe 


of oe 


oo 


pene AXIS ranksbysrank, then the basic solution is 
x = Pi*[ inv(R_11) O J]*Q'*bþ 
[ 0 0 ] 


Ret.: G. H. Golub & C. F. Van Loan, "Matrix Computations", 
Johns Hopkins, 1989. Subsection 5.5.6. 


Per Christian Hansen, UNI-C, 07/11/90. 


x 


Pisce: rank) ~ (Ria): rank, 1: rank) \(O(:;,12rank) '*b) )} 
r 


p - Q*R*Pi*x; 


85 


10. 


Ti; 
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